System and method for the lossless progressive streaming of images over a communication network

ABSTRACT

Disclosed is a lossless image streaming system for the transmission of images over a communication network. The system eliminates the necessity to store a compressed version of the original image, by losslessly streaming ROI data using the original stored image. The imaging system of the present invention also avoids the computationally intensive task of compression of the fall image. When a user wishes to interact with a remote image, the imaging server performs a fast preprocessing step in near real time after which it can respond to any ROI requests also in near real time. When a ROI request arrives at the server, a sophisticated progressive image encoding algorithm is performed, but not for the full image. Instead, the encoding algorithm is performed only for the ROI. Since the size of the ROI is bounded by the size and resolution of the viewing device at the client and not by the size of the image, only a small portion of the full progressive coding computation is performed for a local area of the original image.

[0001] This application claims priority to U.S. Utility patentapplication Ser. No. 09/386,264 filed Aug. 31, 1999, entitled “SYSTEMAND METHOD FOR TRANSMITTING A DIGITAL IMAGE OVER A COMMUNICATIONNETWORK”, and U.S. Provisional Patent Application No. 60/198,017, filedApr. 18, 2000, entitled “LOSSLESS PROGRESSIVE STREAMING OF IMAGES OVERTHE INTERNET”, the entirety of which are both incorporated herein byreference.

BACKGROUND OF THE INVENTION

[0002] 1. Field of the Invention

[0003] This invention relates to systems and methods for transmission ofstill images over relatively low-speed communication channels. Morespecifically the invention relates to progressive image streaming overlow speed communication lines, and may be applied to a variety of fieldsand disciplines, including commercial printing and medical imaging,among others.

BRIEF DESCRIPTION OF THE PRIOR ART

[0004] In a narrow bandwidth environment, a simple transfer to theclient computer of any original image stored in the server's storage isobviously time consuming. In many cases the user only wishes to view alow resolution of the image and perhaps a few more high-resolutiondetails, in these instances it would be inefficient to transfer the fullimage. This problem can be overcome by storing images in some compressedformats. Examples for such formats are standards such as ProgressiveJPEG (W. Pennebaker and J. Mitchel, “JPEG, still image data compressionstandard”, VNR, 1993) or the upcoming JPEG2000 (D. Taubman, “Highperformance scalable image compression with EBCOT”, preprint, 1999).These formats allow progressive transmission of an image such that thequality of the image displayed at the client computer improves duringthe transmission.

[0005] In some application such as medical imaging, it is also necessarythat whenever the user at the client computer is viewing a portion ofthe highest resolution of the image, the progressive streaming willterminate at lossless quality. This means that at the end of progressivetransmission the pixels rendered on the screen are exactly the pixels ofthe original image. The current known “state-of-the-art” waveletalgorithms for progressive lossless streaming all have a major drawback:their rate-distortion behavior is very inferior to the “lossy”algorithms. The implications are serious:

[0006] 1. Whenever the user is viewing any low resolution of the image(at low resolutions the term “lossless” is not well defined) more dataneeds to be sent for the same visual quality.

[0007] 2. During the progressive transmission of the highest resolution,before lossless quality is achieved, more data needs to be sent for thesame visual quality.

[0008] Researchers working in this field are troubled by thesephenomena. As F. Sheng, A. Bilgin, J. Sementilli and M. W. Marcellin sayin [SBSM]: “. . . Improved lossy performance when using integertransforms is a pursuit of our on-going work.” Here is an example: TABLE1 Comparison of the lossy compression performances (implemented by the(7, 9) Wavelet) to a lossless compression (implemented by a reversible(4, 4) Wavelet) of “Barabara”image (PSNR (dB)) ([SBSM]). Rate (bit perpixel) Wavelet 0.1 0.2 0.5 0.7 1.0 Floating Point 7 × 9 24.18 26.6531.64 34.17 36.90 Reversible (4, 4) 23.89 26.41 31.14 33.35 35.65

[0009] As one can see from Table 1, state of the art progressivelossless coding is inferior to lossy coding by more than 1 dB at thehigh bit-rate.

[0010] Indeed, intuitively, the requirement for lossless progressiveimage transmission should not effect the rendering of lower resolutionsor the progressive “lossy” rendering of the highest resolution beforelossless quality is obtained. The final lossless quality should be alayer that in some sense is added to a lossy algorithm with minor (ifany) effect on its performance.

[0011] The main problem with known lossless wavelet algorithms, such asSPIHT [SP1] and CREW [ZASB], is that they use special “Integer ToInteger” transforms (see “Wavelet transforms that map integers tointegers”, A. Calderbank, I. Daubechies, W. Sweldens, B. L. Yeo, J.Fourier Anal. Appl., 1998). These transforms mimic “mathematicallyproven” transforms that work well in lossy compression usingfloating-point arithmetic implementations. But because they areconstraint to be lossless, they do not approximate their floating-pointancestors sufficiently well. Although in all previous work there havebeen attempts to correct this approximation in the progressive codingstage of the algorithm, the bad starting point, an inefficienttransform, prevented previous authors to obtain decent rate-distortionbehavior.

[0012] Our algorithm solves the rate-distortion behavior problem. Usingthe fact that images are two-dimensional signals, we introduce new 2Dlossless Wavelet transforms that approximate much better their lossycounterparts. As an immediate consequence our lossless progressivecoding algorithm has the same rate-distortion of a lossy algorithmduring the lossy part of the progressive transmission.

SUMMARY OF THE INVENTION

[0013] The imaging system that is described below is directed to alossless image streaming system that is different from traditionalcompression systems and overcomes the above problems. By utilizing alossless means of progressive transmission means the pixels rendered onthe screen at the end of transmission are exactly the pixels of theoriginal image that was transmitted. The imaging system disclosed hereineliminates the necessity to store a compressed version of the originalimage, by streaming ROI data using the original stored image. Theimaging system of the present invention also avoids the computationallyintensive task of compression of the full image. Instead, once a userwishes to interact with a remote image, the imaging server performs afast preprocessing step in near real time after which it can respond toany ROI requests also in near real time. When a ROI request arrives atthe server, a sophisticated progressive image-encoding algorithm isperformed, but not for the full image. Instead, the encoding algorithmis performed only for the ROI. Since the size of the ROI is bounded bythe size and resolution of the viewing device at the client and not bythe size of the image, only a small portion of the full progressivecoding computation is performed for a local area of the original image.This local property is also true for the client. The client computerperforms decoding and rendering only for ROI and not the full image.This real time streaming or Pixels-On-Demand™ architecture requiresdifferent approaches even to old ideas. For example, similarly to someprior art, the present imaging system is based on wavelets. But while inother systems wavelet bases are selected according to their codingabilities, the choice of wavelet bases in the present imaging systemdepends more on their ability to perform well in the real timeframework. The system of the present invention supports several modes ofprogressive transmission: by resolution, by accuracy and by spatialorder.

BRIEF DESCRIPTION OF THE DRAWINGS

[0014]FIG. 1 is a system architecture block diagram.

[0015]FIG. 2 is an imaging system workflow diagram.

[0016]FIG. 3 is a flow diagram representing a “lossless progressive byaccuracy” request list for a ROI.

[0017]FIG. 4 is a diagram depicting the client “progressive by accuracy”workflow.

[0018]FIG. 5 is a diagram depicting the server workflow.

[0019]FIG. 6 is a diagram describing the server-preprocessing step.

[0020]FIG. 7 is a diagram describing the low resolution encodingprocess.

[0021]FIG. 8 is a diagram describing the ROI high resolution processing.

[0022]FIG. 9 is a diagram depicting the local forward lossless wavelettransform.

[0023]FIG. 10 is a diagram depicting the local inverse lossless wavelettransform.

[0024]FIG. 11 is a diagram depicting the progressive rendering steps.

[0025]FIG. 12 is a diagram depicting a lossless subband tile wherein thespatial grouping of subband coefficients are at a given resolution andthe halfbit matrix is associated with the hh_(j) subband.

[0026]FIG. 13 a diagram depicting the RGB <->YUV reversible conversion.

[0027]FIG. 14 a diagram depicting a lossless last bit plane data blockin which only hi and hh subband coefficients are scanned.

[0028]FIG. 15 is a sample pseudo-code of the encoding algorithmrepresented by: (a) Least significant bit plane scan pseudo-code (b)Half bit plane scan pseudo-code.

[0029]FIG. 16 is a sample pseudo-code of the decoding algorithmrepresented by: (a) Least significant bit plane scan pseudo-code (b)Half bit plane scan pseudo-code.

[0030]FIG. 17 is a diagram depicting the curve defining the mapping fromOriginal_Image_Depth-bit image to Screen_Depth-bit.

[0031]FIG. 18 is a diagram depicting the decomposition ofone-dimensional signal x to the Low-subband s and the High-subband d andthe separable decomposition of two-dimensional signal X into 4 matrices(subbands): LL,HL,LH and HH.

[0032]FIG. 19 is a diagram depicting the first stage of the 2D separabletransform step in the X-direction.

[0033]FIG. 20 is a diagram depicting the second stage of the 2Dseparable transform step in the Y-direction.

[0034]FIG. 21 is a diagram depicting the application of the full 2DWavelet transform.

[0035]FIG. 22 is a flow diagram representing the least significant bitplane scan of the encoding algorithm.

[0036]FIG. 23 is a flow diagram representing the least significant bitplane scan of the decoding algorithm.

[0037]FIG. 24 is a flow diagram describing a bit plane significance scanof the server-encoding algorithm.

[0038]FIG. 25 is a flow diagram describing the subband tile extractionof the client progressive rendering process.

[0039]FIG. 26 is a diagram depicting the preprocessing multi-resolutionstructure.

DETAILED DESCRIPTION OF THE INVENTION

[0040] 1. Notation and Terminology

[0041] The following notation is used throughout this document. TermDefinition 4D Four dimensional FIR Finite Impulse Response FWT ForwardWavelet Transform GUI Graphical User Interface ID Identification tag IWTInverse Wavelet Transform ROI Region Of Interest URL Uniform ResourceLocator LSB Least Significant Bit RMS Root Square Error FP FloatingPoint PDF Probability Distribution Function

[0042] The following terminology and definitions apply throughout thisdocument. Term Definition Rendering Procedure to output/display a ROI ofan image into a device such as a monitor, printer, etc. MultiresolutionA pyramidal structure representing an image at dyadic resolutions,beginning with the original image as the highest resolution. SubbandTransform/ A method of decomposing an image (time domain) to frequencysubband coefficients components (frequency domain). A representation ofan image as a sum of differences between the dyadic resolutions of theimage's multiresolution. Wavelet Transform/ A special case of SubbandTransform. Wavelet coefficients Progressive Transmitting a given imagein successive steps, where each step Transmission adds more detail tothe rendering of the image Progressive Rendering A sequence of renderingoperations, each adding more detail. Progressive by accuracyTransmit/render strong features first (sharp edges), less significantfeatures (texture) last Progressive by Transmit/render low resolutionfirst, high resolution last resolution Progressive by spatialTransmit/render top of image first, bottom of image last orderDistributed database A database architecture that can be used in anetwork environment Subband/Wavelet tile A group of subband/waveletcoefficients corresponding to a time- frequency localization at somegiven spatial location and resolution/frequency Subband/Wavelet data Anencoded portion of a subband/wavelet tile that corresponds to block someaccuracy layer

[0043] 2. Overview of the Invention

[0044] Referring to FIG. 1, a block diagram is provided depicting thevarious components of the imaging system in one embodiment. A clientcomputer 110 is coupled to a server computer 120 through a communicationnetwork 130.

[0045] In one embodiment, the client computer 110 and server computer120 may comprise a PC-type computer operating with a Pentium-classmicroprocessor, or equivalent. Each of the computers 110 and 120 mayinclude a cache 111, 121 respectively as part of its memory. The servermay include a suitable storage device 122, such as a high-capacity disk,CD-ROM, DVD, or the like.

[0046] The client computer 110 and server computer 120 may be connectedto each other, and to other computers, through a communication network130, which may be the Internet, an Intranet (e.g., a local areanetwork), a wide-area network, or the like. Those having ordinary skillin the art will recognize that any of a variety of communicationnetworks may be used to implement the present invention.

[0047] With reference to FIG. 2, the system workflow is described. Usingany browser type application, the user of the client computer 110connects to the Web Server 140 or directly to the Imaging server 120 asdescribed in FIG. 1. He/she then selects, using common browser tools animage residing on the Image file storage 122. The corresponding URLrequest is received and processed by the Imaging Server 120. In caseresults of previous computations previous computations on the image arenot present in the Imaging Cache 121, the server performs a fastpreprocessing algorithm (see §2.1) in a lossless mode. The result ofthis computation is inserted into the cache 121. Unlike other more“traditional” applications or methods which perform full progressiveencoding of the image in an “offline” type method, the goal of thepreprocessing step is to allow the server, after a relatively fastcomputational step, to serve any ROI specified by the user of the clientcomputer. For example, for a 15M grayscale medical image, using thedescribed (software) server, installed on computer with a Pentiumprocessor, fast disk, running Windows NT, the preprocessing step 501will typically take 3 seconds. This is by order of magnitude faster thana “fill” compression algorithm such as [S], [SP1], [T]. Serving ROI issometimes called “pixels-on-demand” which means a progressivetransmission of any ROI of the image in “real-time”, where the qualityof the view improves with the transfer rate until lossless view isreceived in the client side. Once the preprocessing stage is done, theserver sends to the client a notification that the “image is ready to beserved”. The server also transmits the basic parameters associated withthe image such as dimensions, color space, etc. Upon receiving thisnotification, the client can select any ROI of the image using standardGUI. The ROI is formulated in step 203 into a request list that is sentto the server. Each such request corresponds to a data block (§4). Theorder of requests in the list corresponds to some progressive modeselected in the context of the application such as “progressive byaccuracy” rendering of the ROI. Upon receiving the ROI request list, theserver processes the requests according to their order. For each suchrequest the server checks if the corresponding data block exists in theCache 121. If not, the server then computes the data block, stores it inthe cache and immediately sends it to the client. Once a data block thatwas requested arrives at the client, it is inserted into the cache 111.At various points in time during the transfer process, a decision ruleinvokes a rendering of the ROI. Obviously, if some of the data blocksrequired for a high quality rendering of the ROI, were requested fromthe server, but have not arrived yet, the rendering of the ROI will beof lower quality. But, due to the progressive request order, therendering quality will improve with each received data block in an“optimal” way. In case the user changes the ROI, the rendering task atthe client is canceled and a new request list corresponding to a newROI, sent from the client to the server, will notify the server toterminate the previous computation and transmission task and begin thenew one.

[0048] 3. New Reversible Wavelet Transform

[0049] We will now explain in detail why the rate-distortion behavior ofour progressive lossless algorithm is better than other knownalgorithms. Lossless wavelet transform, must be integer-to-integertransform, such that round-off errors are avoided. In order todemonstrate the difference between lossy and lossless transforms, let uslook at the simplest wavelet, the Haar wavelet. Let x(k) be the k-thcomponent of the one-dimensional discrete signal x. The first forwardHaar transform step, in its accurate “mathematical” form, is defined by:$\begin{matrix}\left\{ \begin{matrix}{{{s(n)} = {\frac{1}{\sqrt{2}}\left( {{x\left( {{2n} + 1} \right)} + {x\left( {2n} \right)}} \right)}},} \\{{{d(n)} = {\frac{1}{\sqrt{2}}\left( {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}} \right)}},}\end{matrix} \right. & (3.1)\end{matrix}$

[0050] Where s is a low-resolution version of x, and d is the“difference” between s and x. In the case of lossless transform,applying the above transform results in round-off error. One possibilityis to apply the transform step suggested by [CDSI]: $\begin{matrix}\left\{ \begin{matrix}{\quad {{{s(n)} = \left\lfloor \frac{{x\left( {{2n} + 1} \right)} + {x\left( {2n} \right)}}{2} \right\rfloor},}} \\{\quad {{d(n)} = {{x\left( {{2n} + 1} \right)} - {{x\left( {2n} \right)}.}}}}\end{matrix} \right. & (3.2)\end{matrix}$

[0051] The symbol └∘┘ is the floor function meaning “greatest integerless than or equal to ∘”, e.g.

└0.5┘=0, └−1.5┘=−1, └2┘=2, └−1┘=−1.

[0052] The one-dimensional transform step is generalized to 2D separabletransform step, by applying the 1D transform step twice, first in theX-direction and than (on the first stage output) in the Y-direction asdescribed in FIG. 18, 19 and 20 (see also [M] 7.7). The full 2D Wavelettransform is applied by using the 2D Wavelet transform step iterativelyin the classic Mallat decomposition of the image (FIG. 21) ([M] 7.7).

[0053] In (3.2) two properties are kept:

[0054] 1. Reversibility, i.e., one can restore x(2n)and x(2n+1), byknowing s(n) and d (n), as follows: $\begin{matrix}\left\{ \begin{matrix}{\quad {{{x\left( {2n} \right)} = {{s(n)} - \left\lfloor \frac{d(n)}{2} \right\rfloor}},}} \\{\quad {{x\left( {{2n} + 1} \right)} = {{d(n)} + {{x\left( {2n} \right)}.}}}}\end{matrix} \right. & (3.3)\end{matrix}$

[0055] 2. De-correlation, i.e., s(n) and d(n) contains the minimalnumber of bits required in order to follow property 1. For example, ifthe transform would have been defined by: $\begin{matrix}\left\{ \begin{matrix}{\quad {{{s(n)} = {{x\left( {{2n} + 1} \right)} + {x\left( {2n} \right)}}},}} \\{\quad {{{d(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},}}\end{matrix} \right. & (3.4)\end{matrix}$

[0056] then the least significant bit of s(n) and d(n) would have beenthe same and saved twice. In other words, there is a correlation betweens(n) and d(n) in (3.4). From the view point of coding this should beavoided since there is a redundancy of transmitting this bit.

[0057] On the other hand, the scaling property, that is a very importantone, is not kept in (3.2). Observe that the value of s(n) computed by(3.2), is smaller than its “real mathematical value” as computed in(3.1), by factor of {square root}{square root over (2)}. Since s(n)should be rounded to an integer number, the fact that s(n) is smallerthan what it should be, increases the round-off error. In lowresolutions, the error is accumulated through the wavelet steps.

[0058] If we take the error as a model of “white noise” added to thei—th resolution in a multi-resolution representation of the image, i.e.X_(i) in FIG. 21, it can be proved that the variance of this noiseexponentially increases as a function of i. This “contamination” to themulti-resolution image damages the coding efficiency at low bit-rates.Let us describe this in detail for the case of the Haar wavelet. We havetwo assumptions in our analysis:

[0059] 1. The parity (least significant bit) of an arbitrary coefficientc, in any of the wavelet steps is a uniformly distributed randomvariable, i.e. $\begin{matrix}{{P\quad {r\left( {c \equiv {0\quad {{mod}2}}} \right)}} = {{P\quad {r\left( {c \equiv {1\quad {{mod}2}}} \right)}} = {\frac{1}{2}.}}} & (3.5)\end{matrix}$

[0060] 2. This parity is independent of other coefficient's parity(Identically independent distributed).

[0061] Our referenced computation, i.e. the accurate computation, is theHaar transform step defined in (3.1). We concentrate on the LL-subbandcoefficients, because the low-resolution subbands are computed fromthem. LL-subband coefficients are the result of a 2D-transform step(FIG. 18)${{{ll}^{({accurate})}\left( {m,n} \right)} = \frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{2m},{{2n} + 1}} \right)} + {x\left( {{{2m} + 1},{2n}} \right)} + {x\left( {{2m},{2n}} \right)}}{2}},.$

[0062] Where m and n are the indices of the row and column of thecoefficient respectively.

[0063] As described in FIG. 18, 19 and 20, according to the step definedin (3.2), we first apply the X-direction step${{s\left( {k,n} \right)} = \left\lfloor \frac{{x\left( {k,{{2n} + 1}} \right)} + {x\left( {k,{2n}} \right)}}{2} \right\rfloor},$

[0064] For each input row x(k,•).

[0065] Under assumption 1 mentioned above, we can write s(k,n) as${{s\left( {k,n} \right)} = {\frac{{x\left( {k,{{2n} + 1}} \right)} + {x\left( {k,{2n}} \right)}}{2} + e}},$

[0066] Where e is a random variable with a probability distributionfunction (P.D.F.) p(•) defined by

[0067] Therefore, $\begin{matrix}\left\{ \begin{matrix}{\quad {{{p\left( {- 0.5} \right)} = {{\Pr \left( {e = {- 0.5}} \right)} = \frac{1}{2}}},}} \\{\quad {{p(0)} = {{\Pr \left( {e = 0} \right)} = {\frac{1}{2}.}}}}\end{matrix} \right. & (3.6) \\{{{E(e)} = {- \frac{1}{4}}},{{{Var}(e)} = {\frac{1}{16}.}}} & (3.7)\end{matrix}$

[0068] We then apply the Y-direction transform by $\begin{matrix}{{{ll}^{({CDSI})}\left( {m,n} \right)} = {\left\lfloor \frac{{s\left( {{{2m} + 1},n} \right)} + {s\left( {{2m},n} \right)}}{2} \right\rfloor = {\frac{{s\left( {{{2m} + 1},n} \right)} + {s\left( {{2m},n} \right)}}{2} + {e^{\prime}.}}}} & (3.8)\end{matrix}$

[0069] As in (3.5) we can represent s(2m+1,n) and s(2m,n) by:$\begin{matrix}{{{s\left( {{{2m} + 1},n} \right)} = {\frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{{2m} + 1},{2n}} \right)}}{2} + e_{1}}},} & (3.9) \\{{s\left( {{2m},n} \right)} = {\frac{{x\left( {{2m},{{2n} + 1}} \right)} + {x\left( {{2m},{2n}} \right)}}{2} + {e_{2}.}}} & (3.10)\end{matrix}$

[0070] Now we can write:

[0071] $\begin{matrix}\begin{matrix}{{l\quad {l^{({C\quad D\quad S\quad I})}\left( {m,n} \right)}} = {\frac{\frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{{2m} + 1},{2n}} \right)}}{2} + e_{1} + \frac{{x\left( {{2m},{{2n} + 1}} \right)} + {x\left( {{2m},{2n}} \right)}}{2} + e_{2}}{2} + e^{\prime}}} \\{= {\frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{{2m} + 1},{2n}} \right)} + {x\left( {{2m},{{2n} + 1}} \right)} + {x\left( {{2m},{2n}} \right)}}{4} + \frac{e_{1}}{2} + \frac{e_{2}}{2} + e^{\prime}}} \\{= {\frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{{2m} + 1},{2n}} \right)} + {x\left( {{2m},{{2n} + 1}} \right)} + {x\left( {{2m},{2n}} \right)}}{4} + e}}\end{matrix} & (3.11)\end{matrix}$

[0072] Where e₁, e₂, e′ are independent (assumption 2 above) randomvariables with${{expectation}\quad \frac{1}{4}},{{Variance}\quad \frac{1}{16}},{{{and}\quad e} = {\frac{e_{1}}{2} + \frac{e_{2}}{2} + {e^{\prime}.}}}$

[0073] Therefore, $\begin{matrix}{{{E(e)} = {{{\frac{1}{2}{E\left( e_{1} \right)}} + {\frac{1}{2}{E\left( e_{2} \right)}} + {E\left( e^{\prime} \right)}} = {- \frac{1}{2}}}},{{{Var}(e)} = {{{Var}\left( {\frac{e_{1}}{2} + \frac{e_{2}}{2} + e^{\prime}} \right)} = {{{\frac{1}{4}{{Var}\left( e_{1} \right)}} + {\frac{1}{4}{{Var}\left( e_{2} \right)}} + {{Var}\left( e^{\prime} \right)}} = {{{\frac{1}{4} \cdot \frac{1}{16}} + {\frac{1}{4} \cdot \frac{1}{16}} + \frac{1}{16}} = {\frac{3}{32}.}}}}}} & (3.12)\end{matrix}$

[0074] Thus, $\begin{matrix}{{l\quad {l^{({C\quad D\quad S\quad I})}\left( {m,n} \right)}} = {\frac{l\quad {l^{({a\quad c\quad c\quad u\quad r\quad a\quad t\quad e})}\left( {m,n} \right)}}{2} + {e.}}} & (3.13)\end{matrix}$

[0075] e represents the approximation error of the LL-subbandcoefficients, results from one 2D transform step. The error relates tothe accurate floating-point computation.

[0076] This was a description of a single 2D-transform step assumingthat the input coefficients are without any error. Now we wish toevaluate the error accumulated after several steps.

[0077] At an arbitrary step i≧0, we can assume that an input coefficientcan be written as:

x _(i)(k,l)=x _(i) ^((accurte))(k,l)+e _(i),

[0078] Where x_(i) ^((accurate))(k,l) is the accurate value achieved byfloating-point computation for all the previous steps, i.e., a stepdefined by $\begin{matrix}\left\{ \begin{matrix}{{{s(n)} = \frac{{x\left( {{2n} + 1} \right)} + {x\left( {2n} \right)}}{2}},} \\{{{d(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},}\end{matrix} \right. & (3.14)\end{matrix}$

[0079] Instead of the integer-to-integer computation in (3.2). Observethat if x_(i) ^((accurate))(k,l) is the i-th resolution imagecoefficient, using (3.14) as the 1D Wavelet step, then $\begin{matrix}{{{x_{i}^{({a\quad c\quad c\quad u\quad r\quad a\quad t\quad e})}\left( {k,l} \right)} = \frac{l\quad {l_{i - 1}^{({a\quad c\quad c\quad u\quad r\quad a\quad t\quad e})}\left( {k,l} \right)}}{2^{i}}},{i \geq 1},} & (3.15)\end{matrix}$

[0080] Where ll_(i−1) ^((accurate))(k,l) is the normalized (L₂−norm)LL-subband coefficient resulting from the i-th 2D transform step using(3.1) as the 1D Wavelet step (see Figure ). e_(i) is the differencebetween x_(i)(k,l) and x_(i) ^((accurate))(k,l) (I.e., the approximationerror of the integer computation made until now). E.g. e₀=0(x₀(k,l) isan original image pixel), while e₁ is a random number with expectation$- \frac{1}{2}$

[0081] And variance $\frac{3}{32}$

[0082] (see (3.12)).

[0083] Using (3.11), we get:

[0084] $\begin{matrix}{{l\quad {l_{i}^{({C\quad D\quad S\quad I})}\left( {m,n} \right)}} = {\frac{{x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + e_{i}^{1} + {x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{2n}} \right)} + e_{i}^{2} + {x_{i}^{({a\quad c\quad c})}\left( {{2m},{{2n} + 1}} \right)} + e_{i}^{3} + {x^{({a\quad c\quad c})}\left( {{2m},{2n}} \right)} + e_{i}^{4}}{4} + e}} \\{= {\frac{{x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{2m},{{2n} + 1}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{2m},{2n}} \right)}}{4} + \frac{e_{i}^{1} + e_{i}^{2} + e_{i}^{3} + e_{i}^{4}}{4} + e}} \\{= {\frac{{x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{2m},{{2n} + 1}} \right)} + {x_{i}^{({a\quad c\quad c})}\left( {{2m},{2n}} \right)}}{4} + e_{i + 1}}} \\{{= {x_{i + 1}^{({a\quad c\quad c\quad u\quad r\quad a\quad t\quad e})} + e_{i + 1}}},}\end{matrix}$

[0085] Where e_(i+1) is defined by${e_{i + 1} = {\frac{e_{i}^{1} + e_{i}^{2} + e_{i}^{3} + e_{i}^{4}}{4} + e}},$

[0086] And corresponds to LL_(i) subband.

[0087] Consequently

E(e _(i+1))=E(e _(i))+E(e),

[0088]${{Var}\left( e_{i + 1} \right)} = {\frac{{Var}\left( e_{i} \right)}{4} + {{{Var}(e)}.}}$

[0089] Observe that${{E(e)} = {- \frac{1}{2}}},{{{Var}(e)} = {\frac{3}{32}.}}$

[0090] As a result, we can write recursive formulas for the errorexpectation and variance after i steps. $\begin{matrix}\begin{matrix}\left\{ \begin{matrix}{\quad {{{E\left( e_{0} \right)} = 0},}} \\{\quad {{{E\left( e_{i + 1} \right)} = {{E\left( e_{i} \right)} - \frac{1}{2}}},}}\end{matrix}\quad \right. \\\left\{ \begin{matrix}{\quad {{{{Var}\left( e_{0} \right)} = 0},}} \\{\quad {{{{Var}\left( e_{i + 1} \right)} = {\frac{{Var}\left( e_{i} \right)}{4} + \frac{3}{32}}},}}\end{matrix} \right.\end{matrix} & (3.16)\end{matrix}$

[0091] The explicit solutions to these formulas are $\begin{matrix}{{{E\left( e_{i} \right)} = {- \frac{i}{2}}},{{{Var}\left( e_{i} \right)} = {\frac{1}{8} - {\frac{1}{2 \cdot 4^{i + 1}}.}}}} & (3.17)\end{matrix}$

[0092] By replacing $\begin{matrix}{{{x_{i}^{({accurate})}\left( {m,n} \right)}\quad {with}\quad \frac{{ll}_{i - 1}^{({accurate})}\left( {m,n} \right)}{2^{i}}\quad {we}\quad {get}}{{{ll}_{i}^{({CDSI})}\left( {m,n} \right)} = {\frac{{ll}_{i}^{({accurate})}\left( {m,n} \right)}{2^{i + 1}} + {e_{i + 1}.}}}} & (3.18)\end{matrix}$

[0093] Thus, the approximation to ll_(i) ^((accurate))(m,n) is

2^(i+1) ll _(i) ^((CDSI))(m,n)=ll _(i) ^((accurate))(m,n)+2^(i+1) e_(i+1).

[0094] The approximation error expectation is${E\left( {2^{i}e_{i}} \right)} = {{2^{i}{E\left( e_{i} \right)}} = {{2^{i}\left( {- \frac{i}{2}} \right)} = {- {{i2}^{i - 1}.}}}}$

[0095] The approximation error variance and standard deviation are${{Var}\left( {2^{i}e_{i}} \right)} = {{4^{i}{{Var}\left( e_{i} \right)}} = {{4^{i}\left( {\frac{1}{8} - \frac{1}{2 \cdot 4^{i + 1}}} \right)} = {{\frac{4^{i} - 1}{8} \approx \frac{4^{i}}{8}} = {2^{{2i} - 3}.{Hence}}}}}$${{Std}\left( {2^{i}e_{i}} \right)} = {\sqrt{{Var}\left( {2^{i}e_{i}} \right)} \approx {\frac{2^{i - 1}}{\sqrt{2}}.}}$

[0096] Let us now evaluate the approximation error of the 3 othersubbands: $\begin{matrix}{{{lh}_{i}^{({CDSI})}\left( {m,n} \right)} = \quad {\left\lfloor \frac{{{ll}_{i - 1}^{({CDSI})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {{ll}_{i - 1}^{({CDSI})}\left( {{{2m} + 1},{2n}} \right)}}{2} \right\rfloor -}} \\{\quad \left\lfloor \frac{{{ll}_{i - 1}^{({CDSI})}\left( {{2m},{{2n} + 1}} \right)} + {{ll}_{i - 1}^{({CDSI})}\left( {{2m},{2n}} \right)}}{2} \right\rfloor} \\{= \quad {\frac{{lh}_{i}^{({accurate})}\left( {m,n} \right)}{2^{i}} + \frac{e_{i}^{1} + e_{i}^{2}}{2} + e^{\prime} - \left( {\frac{e_{i}^{3} + e_{i}^{4}}{2} + e^{''}} \right)}} \\{= \quad {\frac{{lh}_{i}^{({accurate})}\left( {m,n} \right)}{2^{i}} + e_{i}^{L\quad H}}}\end{matrix}$

[0097] Where

[0098] e_(i) ^(k) 1≦k≦4 are identical to the random variable whoseexpectation and variance are given in (3.17).$e_{i}^{L\quad H} = {\frac{e_{i}^{1} + e_{i}^{2}}{2} + e^{\prime} - \left( {\frac{e_{i}^{3} + e_{i}^{4}}{2} + e^{''}} \right)}$

[0099] e′ and e″ are identical to the random variable whose expectationand variance are given in (3.7).

[0100] Thus,

E(e _(i) ^(LH))=0,

[0101] $\begin{matrix}\begin{matrix}{{{Var}\left( e_{i}^{L\quad H} \right)} \approx \quad {{\frac{1}{4}\left( {\frac{1}{8} - \frac{1}{2 \cdot 4^{i + 1}} + \frac{1}{8} - \frac{1}{2 \cdot 4^{i + 1}}} \right)} + \frac{1}{16} +}} \\{\quad {{\frac{1}{4}\left( {\frac{1}{8} - \frac{1}{2 \cdot 4^{i + 1}} + \frac{1}{8} - \frac{1}{2 \cdot 4^{i + 1}}} \right)} + \frac{1}{16}}} \\{= \quad {\frac{1}{4} - {\frac{1}{2 \cdot 4^{i + 1}}.}}}\end{matrix} & (3.19)\end{matrix}$

[0102] The approximation to lh_(i) ^((accurate))(m,n) is

2^(i) lh _(i) ^((CDSI))(m,n)=lh _(i) ^((accurate))(m,n)+2^(i) e _(i)^(LH).

[0103] The approximation error variance and standard deviation are:${{Var}\left( {2^{i}e_{i}^{L\quad H}} \right)} = {{4^{i}{{Var}\left( e_{i}^{L\quad H} \right)}} = {{4^{i} \cdot \left( {\frac{1}{4} - \frac{1}{2 \cdot 4^{i + 1}}} \right)} = {{4^{i - 1} - \frac{1}{8}} \approx {4^{i - 1}.}}}}$

[0104] Therefore

Std(2^(i) e _(i) ^(LH))={square root}{square root over (Var(2^(i) e _(i)^(LH)))}≈{square root}{square root over (4^(i−1))}=2^(i−1).

[0105] A similar approximation error estimation, it can be done with theHL-subband and the HH-subband.

[0106] The approximation error evaluation results are summarized in thefollowing table where the error is the difference between the normalized(in L₂-norm) coefficients according to [CDSI] reversible transform andthe “mathematical” transform (defined in (3.1)). TABLE 2 Normalized (inL₂ − norm) approximation errors of the Wavelet coefficients atresolution i ≧ 0 (i = 0 is the highest resolution) using the (CDSI)reversible Haar transform. Expectation Variance Std LL_(i) − error −(i +1)2^(i) $\frac{4^{i + 1} - 1}{8} \approx 2^{{2i} - 1}$

0.707 · 2^(i) LH_(i) − error 0$\frac{{2 \cdot 4^{i}} - 1}{8} \approx 4^{i - 1}$

0.5 · 2^(i) HL_(i) − error ${- \frac{1}{4}} \cdot 2^{i}$

$\frac{{3 \cdot 4^{i}} - 2}{16} \approx {3 \cdot 4^{i - 2}}$

0.433 · 2^(i) HH_(i) − error 0$\frac{4^{i} - 1}{8} \approx 2^{{2i} - 3}$

0.354 · 2^(i)

[0107] Assuming a low-bit rate transmission where only the coefficientswhose absolute value belongs to the range [2^(b),2^(b+1)) are encoded,for every resolution i, where i is greater than b(less or more). It mustbe noted that the large error implies a significant loss of codingefficiency.

[0108] Instead, we propose a new family of reversible transforms. Theproposed family of integer wavelet transforms has all three properties:

[0109] 1. Reversibility

[0110] 2. De-correlation

[0111] 3. Scaling-i.e. improved approximation of the “mathematical”transform.

[0112] Our 2D transform step is separable also, but the one-dimensionaltransform step, which the 2D transform is based on, is different for theX-direction (step 1901), the Y-direction step applied on the low outputof the X-direction step (step 2001) and the Y-direction step applied onthe high output of the X-direction step (step 2002) as described in FIG.18, 19 and 20.

[0113] The full 2D Wavelet transform is applied by using the 2D Wavelettransform step iteratively in the classic Mallat decomposition of theimage (FIG. 21) ([M] 7.7). As mentioned before, the Wavelet coefficientsin our proposed transform are all scaled, i.e. normalized in L₂-norm asthe Wavelet coefficients computed in the accurate “mathematical”transform.

[0114] In order to achieve the third property (improved approximation ofthe “mathematical” transform), we define an extra matrix we call the“Half bit-matrix” which enables the reversibility of the HighY-transform step (step 2002). The elements that belong to this matrixare bits, such that each bit corresponds to an HH-subband coefficient inthe following interpretation. Let us describe this by the followingexample.

[0115] Supposing

s(n)=7, d ⁽¹⁾(n)=9

[0116] Are a coefficient pair results from a reversible de-correlated1D-wavelet step $\left\{ {\begin{matrix}{{{d^{(1)}(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},} \\{{s(n)} = {\left\lfloor \frac{{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}{2} \right\rfloor.}}\end{matrix}\quad} \right.$

[0117] Now, d⁽¹⁾(n) has to be multiplied by $\frac{1}{2},$

[0118] In order to be scaled.

[0119] The binary form of d⁽¹⁾(n)=9 is

d ⁽¹⁾(n)=1001₂.

[0120] If we now divide d⁽¹⁾(n) by 2 in a floating-point computation weget${d^{F\quad P}(n)} = {{\frac{1}{2}{d^{(1)}(n)}} = {100.1_{2}.}}$

[0121] Let us call the bit, which located on the write-side of thefloating point the “Half Bit”. Observe that the Half Bit of d^(FP)(n) isthe LSB of d⁽¹⁾(n). Therefore, an equivalent way to do this in aninteger computation without loosing the Half-Bit is to calculate firstthe LSB of d⁽¹⁾(n) by

HalfBit(n)=d ⁽¹⁾(n)mod2=9mod2=1,

[0122] Then to shift-write d⁽¹⁾(n) by

d(n)=d ⁽¹⁾(n)>>1=1001>>1=100.

[0123] By saving d(n) and HalfBit(n) we can restore back d⁽¹⁾(n).

[0124] In the proposed transform, this Half-bit is needed in theHH-subband coefficient computation. Therefore in our waveletdecomposition for every HH-subband coefficient (in all scales) there isa corresponding bit, which is the coefficient's Half-bit.

[0125] The Half bit matrix is hidden in the HH-subband in thedescription of FIG. 18, 19 and 20. It is described explicitly in thespecification of the transform as much in the coding algorithm.

[0126] We now present our integer-to-integer versions of the Haartransform and the CDF (1,3) transform for the 2-dimensional case.

[0127] 3.1 Reversible Haar and (CDF) (1,3) Transforms

[0128] 3.1.1Haar Transform

[0129] With respect to FIG. 19:

[0130]3.1.1.1 Step 1901: X-direction

[0131] Forward Step $\begin{matrix}\left\{ \begin{matrix}{{{s(n)} = \left\lfloor \frac{{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}{2} \right\rfloor},} \\{{{d(n)} = {{x\left( {{2n} + 1} \right)} - {{x\left( {2n} \right)}.}}}\quad}\end{matrix} \right. & (3.20)\end{matrix}$

[0132] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{{{{x\left( {2n} \right)} = {{s(n)} - \left\lfloor \frac{d(n)}{2} \right\rfloor}},}\quad} \\{{x\left( {{2n} + 1} \right)} = {{d(n)} + {{x\left( {2n} \right)}.}}}\end{matrix} \right. & (3.21)\end{matrix}$

[0133] With respect to FIG. 20:

[0134] 3.1.1.2 Step 2001: Y-direction—Low Forward Step $\begin{matrix}\left\{ \begin{matrix}{{{{s(n)} = {{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}},}\quad} \\{{{d^{(1)}(n)} = \left\lfloor \frac{{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}{2} \right\rfloor},} \\{{{d(n)} = {2{{d^{(1)}(n)}.}}}\quad}\end{matrix} \right. & (3.22)\end{matrix}$

[0135] Remarks:

[0136] 1. s(n) is a scaled LL-subband coefficient.

[0137] 2. s(n) and d⁽¹⁾(n) are de-correlated and a reversible couple(can be transformed back to x(2n) and x(2n+1)), but d⁽¹⁾(n) is notscaled (it is half its “real value”). Thus, d⁽¹⁾(n) is multiplied by 2.Nevertheless, the LSB of the LH-subband coefficient d(n) is known to be0 and not encoded.

[0138] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{x\left( {{2n} + 1} \right)} = {\frac{1}{2}\left( {{s(n)} + {d(n)} + \left( {{s(n)}{mod}\quad 2} \right)} \right)}},}} \\{\quad {{x\left( {2n} \right)} = {{s(n)} - {{x\left( {{2n} + 1} \right)}.}}}}\end{matrix} \right. & (3.23)\end{matrix}$

[0139] With respect to FIG. 20:

[0140] 3.1.1.3 Step 2002: Y-direction—High Forward Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{d^{(1)}(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},}} \\{\quad {{{{HalfBit}(n)} = {\left( {d^{(1)}(n)} \right){mod}\quad 2}},}} \\{\quad {{{d(n)} - \left\lfloor \frac{d^{(1)}(n)}{2} \right\rfloor},}} \\{\quad {{s(n)} = {{x\left( {2n} \right)} + {{d(n)}.}}}}\end{matrix} \right. & (3.24)\end{matrix}$

[0141] Remark: d⁽¹⁾(n) and s(n) are de-correlated and reversiblecouples, but d⁽¹⁾(n) is not scaled (It is twice its “real value”).Therefore, d⁽¹⁾(n) is divided by 2. By doing that, we lose its leastsignificant bit, which cannot be restored. To solve this problem, asexplained before, we save this bit as the “Half-Bit”. Giving this nameto that coefficient means that its weight is $\frac{1}{2}$

[0142] In the “real mathematical scale”, and it is the least significant(from the approximation point of view).

[0143] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{{{x\left( {2n} \right)} = {{s(n)} - {d(n)}}},} \\{{{d^{(1)}(n)} = {{2{d(n)}} + {{HalfBit}(n)}}},} \\{{x\left( {{2n} + 1} \right)} = {{d^{(1)}(n)} + {{x\left( {2n} \right)}.}}}\end{matrix} \right. & (3.25)\end{matrix}$

[0144] 3.1.2 CDF (1,3) Transform

[0145] 3.1.2.1 Step 1901: X-direction

[0146] With respect to FIG. 19:

[0147] Forward Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{s(n)} = \left\lfloor \frac{{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}{2} \right\rfloor},}} \\{\quad {{{d^{(1)}(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},}} \\{\quad {{d(n)} = {{d^{(1)}(n)} + {\left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{4} \right\rfloor.}}}}\end{matrix} \right. & (3.26)\end{matrix}$

[0148] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{d^{(1)}(n)} = {{d(n)} - \left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{4} \right\rfloor}},}} \\{\quad {{{x\left( {2n} \right)} = {{s(n)} - \left\lfloor \frac{d^{(1)}(n)}{2} \right\rfloor}},}} \\{\quad {{x\left( {{2n} + 1} \right)} = {{x\left( {2n} \right)} + {{d^{(1)}(n)}.}}}}\end{matrix} \right. & (3.27)\end{matrix}$

[0149] With respect to FIG. 20:

[0150] 3.1.2.2 Step 2001: Y-direction—Low Forward Step $\begin{matrix}\left\{ {\begin{matrix}{\quad {{{s(n)} = {{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}},}} \\{\quad {{d^{(1)}(n)} = \left\lfloor \frac{{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)} + \left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{8} \right\rfloor}{2} \right\rbrack}} \\{\quad {{d(n)} = {2{{d^{(1)}(n)}.}}}}\end{matrix},} \right. & (3.28)\end{matrix}$

[0151] Remark: See remarks for (3.22).

[0152] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{s^{(1)}(n)} = {{s(n)} - \left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{8} \right\rfloor}},}} \\{\quad {{{x\left( {{2n} + 1} \right)} = {\frac{1}{2}\left( {{s^{(1)}(n)} + {d(n)} + \left( {{s^{(1)}(n)}{mod}\quad 2} \right)} \right)}},}} \\{\quad {{x\left( {2n} \right)} = {{s(n)} - {{x\left( {{2n} + 1} \right)}.}}}}\end{matrix} \right. & (3.29)\end{matrix}$

[0153] With respect to FIG. 20:

[0154] 3.1.2.3 Step 2002: Y-direction—High Forward Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{s(n)} = \left\lfloor \frac{{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}{2} \right\rfloor},}} \\{\quad {{{d^{(1)}(n)} = {{x\left( {{2n} + 1} \right)} - {x\left( {2n} \right)}}},}} \\{\quad {{d^{(2)}(n)} = {{d^{(1)}(n)} + {\left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{4} \right\rfloor.}}}} \\{{{d(n)} = \left\lfloor \frac{d^{(2)}(n)}{2} \right\rfloor},} \\{{{HalfBit}(n)} = {{d^{(2)}(n)}{mod}\quad 2.}}\end{matrix} \right. & (3.30)\end{matrix}$

[0155] Inverse Step $\begin{matrix}\left\{ \begin{matrix}{\quad {{{d^{(1)}(n)} = {{2{d(n)}} + {{HalfBit}(n)} - \left\lfloor \frac{{s\left( {n - 1} \right)} - {s\left( {n + 1} \right)}}{4} \right\rfloor}},}} \\{\quad {{{x\left( {2n} \right)} = {{s(n)} - \left\lfloor \frac{d^{(1)}(n)}{2} \right\rfloor}},}} \\{\quad {{x\left( {{2n} + 1} \right)} = {{d^{(1)}(n)} + {{x\left( {2n} \right)}.}}}}\end{matrix} \right. & (3.31)\end{matrix}$

[0156] We now compute the approximation error probabilities of ourmethod, and show that it is significantly smaller. We start with theLL-subband error. Assuming e_(i) is the approximation error of theLL-subband in the i-th resolution (FIG. 21), the LL-subband coefficientin the i-th resolution can be written as: $\begin{matrix}\begin{matrix}{{{ll}_{i}^{({new})}\left( {m,n} \right)} = \quad {\left\lfloor \frac{{x_{i}^{({new})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({new})}\left( {{2m},{{2n} + 1}} \right)}}{2} \right\rfloor +}} \\{\quad \left\lfloor \frac{{x_{i}^{({new})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({new})}\left( {{2m},{2n}} \right)}}{2} \right\rfloor} \\{= \quad {\frac{{x_{i}^{({{acc}.})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + e_{i - 1}^{1} + {x_{i}^{({{acc}.})}\left( {{2m},{{2n} + 1}} \right)} + e_{i - 1}^{2}}{2} + e^{\prime} +}} \\{\quad {\frac{{x_{i}^{({{acc}.})}\left( {{{2m} + 1},{2n}} \right)} + e_{i - 1}^{3} + {x_{i}^{({{acc}.})}\left( {{2m},{2n}} \right)} + e_{i - 1}^{4}}{2} + e^{''}}} \\{= \quad {\frac{{x_{i}^{({{acc}.})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({{acc}.})}\left( {{2m},{{2n} + 1}} \right)} + {x_{i}^{({{acc}.})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({{acc}.})}\left( {{2m},{2n}} \right)}}{2} +}} \\{\quad {{\frac{e_{i - 1}^{1} + e_{i - 1}^{2} + e_{i - 1}^{3} + e_{i - 1}^{4}}{2} + e^{\prime} + e^{''}} =}} \\{{= \quad {{{ll}_{i}^{({{acc}.})}\left( {m,n} \right)} + e_{i}}},}\end{matrix} & (3.32)\end{matrix}$

[0157] Where

[0158] ll_(i) ^((new))(m,n) is the new transform LL-subband coefficient(i-th resolution).

[0159] e_(i−1) ^(k) for 1≦k≦4, are identical random variablerepresenting the error from the previous level.

[0160] e′ and e″ are random variables with P.D.F. defined in (3.6).

[0161] X_(i) ^((acc.))(m,n) is the i-th resolution image coefficientusing (3.1) as the 1D Wavelet step.

[0162] ll_(i) ^((acc.))(m,n) is the normalized (L₂-norm) LL-subbandcoefficient resulting from the i-th 2D transform step using (3.1) as the1D Wavelet step (see FIG. 21).$e_{i} = {\frac{e_{i - 1}^{1} + e_{i - 1}^{2} + e_{i - 1}^{3} + e_{i - 1}^{4}}{2} + e^{\prime} + {e^{''}.}}$

[0163] Consequently $\begin{matrix}{{{E\left( e_{i} \right)} = {{\frac{4{E\left( e_{i - 1} \right)}}{2} + \left( {- \frac{1}{4}} \right) + \left( {- \frac{1}{4}} \right)} = {{2{E\left( e_{i - 1} \right)}} - \frac{1}{2}}}},} & (3.33) \\{{{Var}\left( e_{i} \right)} = {{\frac{4{{Var}\left( e_{i - 1} \right)}}{4} + \frac{1}{16} + \frac{1}{16}} = {{{Var}\left( e_{i - 1} \right)} + {\frac{1}{8}.}}}} & (3.34)\end{matrix}$

[0164] By knowing that e⁻¹=0 we get $\begin{matrix}\left\{ \begin{matrix}{{{E\left( e_{i} \right)} = {\frac{1}{2} - 2^{i}}},} \\{{{Var}\left( e_{i} \right)} = {\frac{i + 1}{8}.}}\end{matrix} \right. & (3.35)\end{matrix}$

[0165] Now we can easily evaluate the approximation error of the 3 othersubbands:

lh_(i) ^((new))(m,n) $\begin{matrix}\begin{matrix}{{{lh}_{i}^{({new})}\left( {m,n} \right)} = \quad {2\left\lfloor \frac{\begin{matrix}{\left\lfloor \frac{{x_{i}^{({new})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({new})}\left( {{2m},{{2n} + 1}} \right)}}{2} \right\rfloor -} \\\left\lfloor \frac{{x_{i}^{({new})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({new})}\left( {{2m},{2n}} \right)}}{2} \right\rfloor\end{matrix}}{2} \right\rfloor}} \\{= \quad {\frac{\begin{matrix}{\left( {{x_{i}^{({{acc}.})}\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x_{i}^{({{acc}.})}\left( {{2m},{{2n} + 1}} \right)}} \right) -} \\\left( {{x_{i}^{({{acc}.})}\left( {{{2m} + 1},{2n}} \right)} + {x_{i}^{({{acc}.})}\left( {{2m},{2n}} \right)}} \right)\end{matrix}}{2} +}} \\{\quad {\frac{e_{i - 1}^{1} + e_{i - 1}^{2} - e_{i - 1}^{3} - e_{i - 1}^{4}}{2} + e^{\prime} - e^{''} + {2e^{\prime\prime\prime}}}} \\{{= \quad {{{lh}_{i}^{({{acc}.})}\left( {m,n} \right)} + e_{i}^{LH}}},}\end{matrix} & (3.36)\end{matrix}$

 =lh _(i) ^((acc.))(m,n)+e _(i) ^(LH),

[0166] Where

[0167] lh_(i) ^((new))(m,n) is the new transform LH-subband coefficient(i-th resolution).

[0168] lh_(i) ^((accurate))(m,n) is the normalized (L₂-norm) LH-subbandcoefficient resulting from the i-th 2D transform step using (3.1) as the1D Wavelet step (see FIG. 21).

[0169] e′″ is a random variable with P.D.F. defined in (3.6).$e_{i}^{LH} = {\frac{e_{i - 1}^{1} + e_{i - 1}^{2} - e_{i - 1}^{3} - e_{i - 1}^{4}}{2} + e^{\prime} - e^{''} + {2e^{\prime\prime\prime}}}$

[0170] all other symbols are defined like in (3.32).

[0171] Hence $\begin{matrix}{{{E\left( e_{i}^{LH} \right)} = {{\frac{{2{E\left( e_{i - 1} \right)}} - {2{E\left( e_{i - 1} \right)}}}{2} + \left( {- \frac{1}{4}} \right) - \left( {- \frac{1}{4}} \right) + {2 \cdot \left( {- \frac{1}{4}} \right)}} = {- \frac{1}{2}}}},} & (3.37) \\{{{Var}\left( e_{i}^{LH} \right)} = {{\frac{4{{Var}\left( e_{i - 1} \right)}}{4} + \frac{1}{16} + \frac{1}{16} + {4 \cdot \frac{1}{16}}} = {\frac{i + 3}{8}.}}} & (3.38)\end{matrix}$

[0172] Similar estimation can be done for the HL and the HH subbands.

[0173] The error estimation (for all subbands) are summarized in thefollowing table where the error is the difference between the normalized(in L₂-norm) coefficients according to our new reversible transform andthe “mathematical” transform (defined in (3.1)). TABLE 3 Normalized (inL₂ − norm) approximation errors of the Wavelet coefficients atresolution i ≧ 0 (i = 0 is the highest resolution) using the proposedreversible Haar transform. The result for the LL-subband is valid forthe proposed reversible (1,3) transform also. Expecation Variance StdLL_(i) − error $\frac{1}{2} - 2^{i}$

$\frac{i + 1}{8}$

$\sqrt{\frac{i + 1}{8}}$

LH_(i) − error $- \frac{1}{2}$

$\frac{i + 3}{8}$

$\sqrt{\frac{i + 3}{8}}$

HL_(i) − error $- \frac{1}{4}$

$\frac{i}{8} + \frac{1}{16}$

$\sqrt{\frac{i}{8} + \frac{1}{16}}$

HH_(i) − error $- \frac{1}{4}$

$\frac{i}{8} + \frac{1}{16}$

$\sqrt{\frac{i}{8} + \frac{1}{16}}$

[0174] The meaning of this result is that in a low bit-rate, where onlylarge coefficients are encoded, this error is negligible.

[0175] 4. Imaging Protocol and Distributed Database

[0176] Dividing the data into tiles and bit-planes

[0177] For the purpose of efficient rendering the coefficients may besub-divided into tiles. The tiles of this invention differ from previousart as shown in FIG. 12. As in the lossy algorithm, here also thesubband tiles are further decomposed to subband data blocks. Each datablock of lossless subband tile (FIG. 12) will have a 4D coordinate

(t_x, t_y, t_resolution,t_bitPlane)

[0178] Where 0≦t_bitPlane≦maxBitPlane(t_resolution).

[0179] Each data block contains the following data in encoded format:

[0180] 1. For t_(—bitPlane≧)2:

[0181] a. A list of the indices of all subband coefficients whoseabsolute value is in the range [2^(t) ^(_(—)) ^(bitPlane−1), 2^(t)^(_(—)) ^(bitPlane)).

[0182] b. The sign of all the coefficients of a.

[0183] c. For t_bitPlane>2, an additional precision bit for anycoefficient that belongs to the current bit plane or any higher bitplane.

[0184] 2. For t_bitPlane=1, which we call the “least significant bitplane”:

[0185] a. A list of the indices of HL-subband and HH-subbandcoefficients whose absolute value belongs to the set {−1,0,1}.

[0186] b. A “zero flag” for each coefficient of a, which indicates ifthe coefficient is equal to zero or not.

[0187] c. The sign of all the coefficients of a, whose “zero flag” isfalse.

[0188] d. The LSB of the HL-subband and HH-subband coefficients thatbelong to higher bit plane.

[0189] Remark:

[0190] Since the LH-subband contains only even coefficients, their LSBmust be zero and is not coded.

[0191] 3. For t_bitPlane=0, which we call the “half bit plane”, a matrixof $\left( \frac{tileLength}{2} \right)^{2}$

[0192] Bits associated with HH-subband coefficients as their last “halfbit” (See (3.24) or (3.30)).

[0193] 5. The Progressive Subband Coding Algorithm

[0194] 5.1 The Encoding Algorithm

[0195] The encoding algorithm of the present invention is performed atthe server 120. In the present imaging system this rather time consumingtask is performed locally in near real-time for a ROI, and not on thefull image. The encoding algorithm is described for images with a singlecolor component, such as grayscale images, but of course may also beapplied to images with multiple color components. The straightforwardgeneralization for an arbitrary number of components will be explainedlater.

[0196] The lossless algorithm receive as input the following parameters:TABLE 4 Lossless Encoding Algorithm Input Parameters Variable Meaningcoef Matrix of subband coefficients, containing$3 \times \left( \frac{tileLength}{2} \right)^{2}{coefficients}$

HalfBit Matrix of bits containing$\left( \frac{tileLength}{2} \right)^{2}$

bits.

[0197] The coding strategy is similar in some sense to that described inA. Said and W. Pearlman, “A new, fast and efficient image codec based onset partitioning”, IEEE Trans. Circuits and Systems for video Tech.,Vol. 6, No. 3, pp. 243-250, 1996, but the preferred embodiment uses no“Zero Tree” data. For all the data blocks with t_bitplane ≧2, we use thelossy encoding algorithm described in previous art with the parameters:

[0198] coef:=coef (The lossy parameter coef initialized with thelossless parameter coef)

[0199] equalBinSize:=True

[0200] ε_(c):=2

[0201] Remark: The lossy algorithm encodes all the bit-plane informationfor t_bitplane≧2.

[0202] For t_bitPlane≦1, i.e. the least significant bit plane (of thelossless algorithm) and the half bit plane, we use a different algorithmdescribed in 5.1.3.

[0203] 5.1.1 Encoding Algorithm Initialization

[0204] The lossless encoding algorithm initialization is the same as thelossy algorithm of § 4.1.1 in the above-cited Ser. No. 09/386,264, whichdisclosure is incorporated herein by reference. In order to initializethe encoding algorithm, the following procedure is performed:

[0205] 1. Assign to each coefficient coef (x,y) its bit plane b(x,y)such that:

|coef(x,y)|∈[ε_(c)2^(b),ε_(c)2^(b+1))

[0206] 2. Compute the maximum bit plane over all such coefficients:${{maxBitPlane}({tile})} = {\max\limits_{x,y}\left( {b\left( {x,y} \right)} \right)}$

[0207] 3. Write the value of maxBitPlane(tile) using one byte as theheader of the data block:

(t_x,t_y,t_resolution, maxBitPlane(t_resolution))

[0208] 4. Initialize all the coefficients as members of theircorresponding Type16 group.

[0209] 5. Initialize a list of significant coefficients to be empty.

[0210] 6. Initialize a coefficient approximation matrix coef as zero.

[0211] 5.1.2 The Outer Loop

[0212] The outer loop of the encoding algorithm scans the bit planesfrom b=maxBitPlane(tile) to b=0. The output of each such bit plane scanis the subband data block. Since the last stage of the encodingalgorithm is arithmetic encoding of given symbols, at the beginning ofeach scan the arithmetic encoding output module is redirected to thestorage area allocated for the particular data block. Once the bit planescan is finished and the data block has been encoded, the output streamis closed and the bit plane b is decremented. After the outer loop isfinished the following stages are performed:

[0213] 1. Least significant bit plane is encoded (t_bitPlane=1).

[0214] 2. Half bit plane is encoded (t_bitPlane=0).

[0215] The output of the least significant bit plane scan is the datablock (FIG. 14):

(t_x, t_y, t_resolution, t_bitPlane=1).

[0216] The half bit plane data block is:

(t_x, t_y, t_resolution, t_bitPlane=0).

[0217] 5.1.3 Bit-plane Scan

[0218] For t_bitPlane≧2, the framework of the bit plane scan isdescribed in FIG. 24, while the pseudo code is given in the above-citedSer. No. 09/386,264, which disclosure is incorporated herein byreference. The scan, for a given level b (b≧2), encodes all of thecoefficients'data corresponding to the absolute value interval[2^(b−1),2^(b)).

[0219] Remark: The encoder method isLastBitPlane( ) is associated to thet_bitPlane=2. For the least significant bit plane, a pseudo code isdescribed in FIG. 15, while a flow chart is described in FIG. 22.

[0220] Explanations to the least significant bit encoding algorithm:

[0221] 1. The coefficients scanning procedure, i.e. moreCoef( )procedure in FIG. 15 or “More coefficients?” in FIG. 22 includes all thecoefficients belong to the HH and the HL-subband (FIG. 14). TheLH-subband is skipped, since the least significant bit of eachcoefficient in it is zero (see remark 2 for (3.22)).

[0222] 2. The procedure isCoefReporeted( ) (“Is coefficient reported?”in the flow chart) returns false if the coefficient is one of {−1, 0,1},i.e. in all higher bit plane scans it was insignificant, otherwise itreturns true.

[0223] 3. The procedure isCoefExactZero( ) (“Coefficient is zero?” inthe flow chart) returns true iff the coefficient is zero.

[0224] 4. The procedure getCoefSign( ) returns the coefficient's sign.

[0225] For the half bit plane, a pseudo code is described in FIG. 16.

[0226] 5.2 The Decoding Algorithm

[0227] Obviously, this algorithm is a reversed step of the encodingalgorithm of section 5.1, performed in the server 120. The clientcomputer 110 during the progressive rendering operation performs thedecoding algorithm. Similar to the encoding algorithm, the decodingalgorithm is described for an image with one component (such as agrayscale image), but of course could also be used with an image withmore than one component. The input parameters to the lossless algorithmare given below: TABLE 5 Lossless decoding algorithm input parametersVariable Meaning coef Empty matrix of subband coefficients to be filledby the decoding algorithm. HalfBit Matrix of bits containing$\left( \frac{tileLength}{2} \right)^{2}$

bits.

[0228] For all the data blocks with t_bitplane≧2, a “lossy” decodingalgorithm is utilized. The input parameters for the lossy algorithm are:

[0229] coef:=coef

[0230] equalBinSize:=True

[0231] ε_(c):=2

[0232] 5.2.1 Decoding Algorithm Initialization

[0233] 1. Assign to each coefficient Coef (x,y) the value zero.

[0234] 2. Assign to each bit belongs to the HalfBit matrix the valuezero.

[0235] 3. Initialize all the coefficients as members of theircorresponding Type16 group.

[0236] 4. Initialize the list of significant coefficients to be empty.

[0237] 5. If the “first” data block

(t_x,t_y,t_resolution, maxBitPlane(t_resolution))

[0238] Is available at the client, read the first byte, which is thevalue of maxBitPlane (tile).

[0239] 5.2.2 The Outer Loop

[0240] Upon the completion of the outer loop in 5.1.2 the followingstages are preformed:

[0241] 1. The decoding algorithm scans the least significant bit plane.The input to this stage is encoded data block (t_x, t_y, t_resolution,LeastSiginificant_bitPlane).

[0242] 2. The decoding algorithm scans the half bit plane. The input tothis stage is encoded data block (t_x, t_y, t_resolution,Half_bitPlane).

[0243] 5.2.3 Bit Plane Scan

[0244] The preferred embodiment follows the lossy prior art of theabove-cited Ser. No. 09/386,264, which disclosure is incorporated hereinby reference, for t_bitPlane≧2. The scan, for a given level b, decodesall of the coefficients' data corresponding to the absolute valueinterval [ε2^(b),ε2^(b+1)). Pseudo codes of the least significant bitplane scan and half bit plane scan are described in FIG. 16.

[0245] Explanations to the pseudo code:

[0246] 1. The decoder's method morecoef( ) scans all the coefficients inthe HH, HL and LH subband. But, since the LH-subband is skipped in theencoding algorithm, we don't call to decodeSymbol( ) for itscoefficients. Instead of this, we update their least significant bit aszero.

[0247] 2. Recall that LH-subband coefficients that have not beenreported until the least significant bit-plane must be zero since theyare known to be even.

[0248] 6. Client Workflow

[0249] With reference to FIG. 4, we describe the workflow at the clientunit 110. Any new ROI generated by the user's action such as a zoom-in,a scroll, or a luminance tuning invokes in step 401 a call from the GUIinterface to the client imaging module with the new ROI view parameters.The client imaging module then computes in step 402 which data blocksare required for the rendering of the ROI and checks if these datablocks are available in the client cache. If not, their coordinate isadded to a request list ordered according to some progressive mode. Therequest list is then encoded in step 403 and sent to the server. Theserver responds to the request by sending back to the client a stream ofdata blocks, in the order in which they were requested. In step 404 theclient inserts them to their appropriate location in the distributeddatabase. At various points in time step 405, a rendering of the ROI, isinvoked. Naturally, the rendering operation is progressive and is ableto use only the currently available data in the client's database.

[0250] 6.1 Step 401: Receiving the ROI Parameters

[0251] The imaging module on the client computer 120 receives from theGUI interface module view parameters detailed in Table 5. Theseparameters are used to generate a request list for the ROI. The sameparameters are used to render the ROI. TABLE 5 Variable MeaningworldPolygon A 2D polygon in the original image coordinate system scaleThe view resolution. 0 < scale < 1 implies a low view resolution, scale= 1 original resolution and scale > 1 a higher than original viewresolution deviceDepth A number in the set {8, 16, 24} representing thedepth of the output device (screen, printer) viewQuality A qualityfactor in the range [1,7] where 1 implies very low quality and 7 implieslossless quality luminanceMap If active: a curve defining a mapping ofmedical images with more than 8 bits (typically 10, 12, 16 bits) pergrayscale values to an 8 bit screen progressiveMo One of: Progressive ByAccuracy, Progressive By Resolution, Progressive by Spatial Order

[0252] The basic parameters of a ROI are worldPolygon and scale whichdetermine uniquely the ROI view. If the ROI is to be rendered onto aviewing device with limited resolution, then a worldPolygon containing alarge portion of the image will be coupled by a small scale. In the casewhere the rendering is done by a printer, the ROI could be a strip of aproof resolution of the original image that has arrived from the servercomputer 120. This strip is rendered in parallel to the transmission,such that the printing process will terminate with the end oftransmission. The other view parameters determine the way in which theview will be rendered. The parameters deviceDepth and viewQualitydetermine the quality of the rendering operation. In cases the viewingdevice is of low resolution or the user sets the quality parameter to alower quality, the transfer size can be reduced significantly.

[0253] The parameter luminanceMap is typically used in medical imagingfor grayscale images that are of higher resolution than the viewingdevice. Typically, screens display grayscale images using 8 bits, whilemedical images sometimes represent each pixel using 16 bits. Thus, it isnecessary to map the bigger pixel range to the smaller range of [0,255].

[0254] Lastly, the parameter progressiveMode determines the order inwhich data blocks should be transmitted from the server 120. The“Progressive By Accuracy” mode is the best mode for viewing in lowbandwidth environments. “Progressive By Resolution” mode is easier toimplement since it does not require the more sophisticated accuracy (bitplane) management and therefore is commonly found in other systems. Thesuperiority of the “progressive by accuracy” mode can be mathematicallyproven by showing the superiority of “non-linear approximation” over“linear approximation” for the class of real-life images. See, e.g., R.A. DeVore, “Nonlinear approximation”, Acta Numerica, pp. 51-150, 1998.

[0255] The “Progressive by Spatial Order” mode is designed, for example,for a “print on demand” feature where the ROI is actually a lowresolution “proof print” of a high resolution graphic art work. In thismode the image data is ordered and received in a top to bottom order,such that printing can be done in parallel to the transmission.

[0256] However, since lossless compression is mostly required in medicalimages transmission, where typically more than 8 bits images are used,we amplify our discussion on the curve (luminanceMap hereinabove) whichdefines the mapping from the original image gray scale range (typically10,12,16 bits) to an 8-bit screen. Further more, in medical imagesviewing, regardless of the original image depth, mapping is required inorder to control the brightness and contrast of the image.

[0257] 6.1.1 Luminance Mapping

[0258] Mapping from original image depth (e.g. 10,12,16 bits ) to screendepth (typically 8-bits), is defined by a monotonic function (FIG. 17):

ƒ:[0,2^(original) ^(_(—)) ^(image) ^(_(—)) ^(depth)−1]→[0,2^(screen)^(_(—)) ^(depth)−1].  (6.1)

[0259] The curve influences not only the mapping, i.e. the drawing tothe screen, but also the request from the server. To understand that,let us concentrate in the maximal gradient of the curve (FIG. 17). In alossy mode, the request is created such that the image approximation inthe client side is close enough to the original image, i.e., the RMS(Root Mean Square Error) is visually negligible. When a curve (mappingfunction) is applied, the RMS can be increased or reduced. The maximalRMS increasing factor depends on the maximal gradient of the curve asfollows: $\begin{matrix}{{{{RMS\_ increasing}{\_ factor}} = {\frac{{RMS}\left( {{f(I)},{f\left( \hat{I} \right)}} \right)}{{RMS}\left( {I,\hat{I}} \right)} \leq {\max \left( f^{\prime} \right)}}},} & (6.2)\end{matrix}$

[0260] Where

[0261] I is the original image

[0262] Î is the approximated image

[0263] ƒ is the mapping function${{RMS}\left( {I_{1},I_{2}} \right)} = \frac{{{I_{1} - I_{2}}}_{L_{2}}}{Image\_ size}$

[0264] max (ƒ′) is the maximal gradient of the curve.

[0265] We consider the worst case of the RMS increasing factor i.e.:

RMS_increasing_factor=Maximal_gradient=max(ƒ′)

[0266] If the RMS increasing factor is greater than 1, it means that the“new RMS” may be greater than we consider as visually negligible error.Thus, the request list should be increase (more bit-planes should berequested from the server) in order to improve the approximationaccuracy. Conversely, if the RMS increasing factor is smaller than 1,the request listing can be reduced. The exact specification of this isgiven in the following section.

[0267] 6.2 Step 402: Creating the Request List

[0268] In step 402 using the ROI view parameters, the client imagingmodule at the client computer 110 calculates data blocks request listordered according to the progressiveMode. Given the parametersworldPolygon and Scale, it may be determined which subband tiles in the“frequency domain” participate in the reconstruction of the ROI in the“time domain”. These tiles contain all the coefficients that arerequired for an “Inverse Subband/Wavelet Transform” (IWT) step thatproduces the ROI. First, the parameter dyadicResolution (ROI) iscomputed, which is the lowest possible dyadic resolution higher than theresolution of the ROI. Any subband tiles of a higher resolution thandyadicResolution (ROI) do not participate in the rendering operation.Their associated data blocks are therefore not requested, since they arevisually insignificant for the rendering of the ROI. If scale≧1, thenthe highest resolution subband tiles are required. Ifscale≦2^(1−numberOfResolutions) then only the lowest resolution tile isrequired. For any other value of scale we perform the mapping describedbelow in Table 6. TABLE 6 scale highestSubbandResolution scale ≦2^(1 − numberOfResolutions) 1 2^(1 −numberOfResolutions) > scale ≦1numberOfResolutions − └− log₂ (scale)┘ scale > 1 numberOfResolutions

[0269] Once it has been determined which subband tiles participate inthe rendering of the ROI, it is necessary to find which of their datablocks are visually significant and in what order they should berequested. Using well known rate/distortion rules from the field ofimage coding (such as is described in S. Mallat and F. Falzon,“Understanding image transform codes”, Proc. SPIE Aerospace Conf.,1997), it is not too difficult to determine an optimal order in whichthe data blocks should be ordered by the client imaging module (and thusdelivered by the server 120). This optimal order is described in steps301-310 of FIG. 3 for the “Progressive By Accuracy” mode. The underlyingmathematical principal behind this approach is “Non-LinearApproximation”.

[0270] First, the subband coefficients with largest absolute values arerequested since they represent the most visually significant data suchas strong edges in the image. Notice that high resolution coefficientswith large absolute value are requested before low resolutioncoefficients with smaller absolute value. Within each given layer ofprecision (bit plane) the order of request is according to resolution;low resolution coefficients are requested first and the coefficients ofhighestSubbandResolution are requested last.

[0271] The main difficulty of this step is this: Assume a subband tileis required for the rendering of the ROI. This means thatt_resolution≦dyadicResolution (ROI) and the tile is required in the IWTprocedure that reconstructs the ROI. It must be understood which of thedata blocks associated with the subband tile represent visuallyinsignificant data and thus should not be requested. Sending all of theassociated data blocks will not affect the quality of the progressiverendering. However, in many cases transmitting the “tail” of data blocksassociated with high precision is unnecessary since it will be visuallyinsignificant. In such a case, the user will see that the transmissionof the ROI from the server 120 is still in progress, yet the progressiverendering of the ROI seems to no longer to change the displayed image.

[0272] Additionally, the influence of the luminance mapping on theaccuracy level of the requested data block as described below. Supposingfor some t_x,t_y and t_resolution, the set

{(t_x,t_y,t_resolution,t_bitplane)|T≦t_bitPlane≦maxPlaneBit(t_resolution)}

[0273] Is requested where T is the minimal bit plane required to thecurrent view. Here, where the luminance mapping is taken in account, thevalue of T might be increased or decreased.

[0274] The number of bit planes reduced (added) from the request list is$\left\lfloor {\log_{2}\left( \frac{1}{Maximal\_ gradient} \right)} \right\rfloor.$

[0275] i.e., for those t_x,t_y and t_resolution mentioned before, thefollowing set is requested:

{(t_x,t_y,t_resolution, t_bitPlane)|T′≦t_bitPlane≦maxPlaneBit(t_resolution)}

[0276] Where$T^{\prime} = {T + {\left\lfloor {\log_{2}\left( \frac{1}{Maximal\_ gradient} \right)} \right\rfloor.}}$

[0277] Examples:

[0278] 1. Given

[0279] Image depth of 12-bit

[0280] Screen depth of 8-bit

[0281] Linear luminance mapping, I.e.,${{Maximal}\quad {gradient}} = {\frac{2^{8}}{2^{12}} = {2^{- 4}.}}$

[0282] The number of bit planes reduced from the request list is:$\left\lfloor {\log_{2}\left( \frac{1}{{Maximal}\quad {gradient}} \right)} \right\rfloor = {\left\lfloor {\log_{2}\left( \frac{1}{2^{- 4}} \right)} \right\rfloor = 4.}$

[0283] 2. Given a luminance mapping with Maximal gradient=2

[0284] The number of bit planes reduced from the request list is:$\left\lfloor {\log_{2}\left( \frac{1}{{Maximal}\quad {gradient}} \right)} \right\rfloor = {\left\lfloor {\log_{2}\left( \frac{1}{2} \right)} \right\rfloor = {- 1.}}$

[0285] Thus one bit plane is added to the original set.

[0286] 6.3 Step 403: Encoding the Request List

[0287] The client imaging module in the client computer 110 encodes therequest list into a request stream that is sent to the server computer120 via the communication network 130 (FIG. 1). The request listencoding algorithm is a simple rectangle-based procedure. The heuristicsof the algorithm is that the requested data block usually can be groupedinto data block rectangles. From the request list of data blocks indexedby the encoding algorithm computes structures of the type

{(t_x,t_y,t_resolution, t_bitPlane),n_(x),n_(y}, n) _(x),n_(y)≧1  (1.3)

[0288] Each such structure represents the n_(x)xn_(y) data blocks

{(t_x+i,t_y+j,t_resolution,t_bitplane)}, 0≦i<n_(x), 0≦j<n_(y)

[0289] The encoding algorithm attempts to create the shortest possiblelist of structures, collecting the data blocks to the largest possiblerectangles can do this. It is important to note that the algorithmensures the order of data blocks in the request list is not changed,since the server 120 will respond to the request stream by transmittingdata blocks in the order in which they were requested. A good example ofwhen this works well is when a user zooms in into a ROI at a highresolution that was never viewed before. In such a case the request listmight be composed of hundreds of requested data blocks, but they will becollected to one (x, y) rectangle for each pair (t_resolution,t_bitPlane).

[0290] 6.4 Step 404: Receiving the Data Blocks

[0291] The client computer 110 upon receiving from the server computer120 an encoded stream containing data blocks, decodes the stream andinserts the data blocks into their appropriate location in thedistributed database using their ID as a key. The simple decodingalgorithm performed here is a reversed step of the encoding infra. Sincethe client 110 is aware of the order of the data blocks in the encodedstream, only the size of each data block need be reported along with theactual data. In case the server 120 informs of an empty data block, thereceiving module marks the appropriate slot in the database as existingbut empty.

[0292] Recall that the subband tile associated with each data block isdenoted by the first three coordinates of the four coordinates of a datablock (t_x,t_y,t_resolution). From the subband tile's coordinates thedimensions are calculated of the area of visual significance; that is,the portion of the ROI that is affected by the subband tile. Assume thateach subband tile is of length tileLength and that the wavelet basisused has a maximal filter size maxFilterSize, then defining hFilterSize:=┌maxFilterSize/2┐ and factor:=numberOfResolutions—t_resolution+1, wehave that the dimensions of the affected region of the ROI (in theoriginal image's coordinate system) are

[t_x×tileLength^(factor)-hFilterSize^(factor),(t_x+1)×tileLength^(factor)+hFilterSize^(factor)]×

[t_y×tileLength^(factor)-hFilterSize^(factor),(t_y+1)×tileLength^(factor)+hFilterSize^(factor])

[0293] These dimensions are merged into the next rendering operation'sregion. The rendering region is used to efficiently render only theupdated portion of the ROI.

[0294] 6.5 Progressive Rendering

[0295] During the transmission of ROI data from the server to theclient, the client performs rendering operations of the ROI. To ensurethat these rendering tasks do not interrupt the transfer, the clientruns two program threads: communications and rendering. The renderingthread runs in the background and draws into a pre-allocated“off-screen” buffer. Only then does the client use device and systemdependant tools to output the visual information from the “off-screen”to the rendering device such as the screen or printer.

[0296] The rendering algorithm performs reconstruction of the ROI at thehighest possible quality based on the available data at the client. Thatis, data that was previously cached or data that “just” arrived from theserver. For efficiency, the progressive rendering is performed only forthe portion of the ROI that is affected by newly arrived data.Specifically, data that arrived after the previous rendering task began.This “updated region” is obtained using the method of step 404 describedin §6.4.

[0297] The parameters of the rendering algorithm are composed of twosets:

[0298] 1. The ROI parameters described in Table 3.

[0299] 2. The parameters transmitted from the server explained in Table5, with the exception of the jumpsize parameter, which is a “serveronly” parameter.

[0300] The rendering algorithm computes pixels at the dyadic resolutiondyadicRresolution(ROI). Recall that this is the lowest possible dyadicresolution that is higher than the resolution of the ROI. The obtainedimage is then resized to the correct resolution. Using a tiling of themultiresolution representation of the ROI, the steps of the algorithmare performed on a tile by tile basis as described in FIG. 11. Since thetiles' length are tileLength, which is typically chosen as 64, therendering algorithm is memory efficient.

[0301] 6.5.1 The Rendering Rate

[0302] As ROI data is transmitted to the client 110, the renderingalgorithm is performed at certain time intervals of a few seconds. Ateach point in time, only one rendering task is performed for any givendisplayed image. To ensure that progressive rendering does not become abottleneck, two rates are measured: the data block transfer rate and theROI rendering speed. If it predicted that the transfer will be finishedbefore a rendering task, a small delay is inserted, such that renderingwill be performed after all the data arrives. Therefore, in a slownetwork scenario (as the Internet often is), for almost the entireprogressive rendering tasks, no delay is inserted. With the arrival ofevery few kilobytes of data, containing the information of a few datablocks, a rendering task visualizes the ROI at the best possiblequality. In such a case the user is aware that the bottleneck of the ROIrendering is the slow network and has the option to accept the currentrendering as a good enough approximation of the image and not wait forall the data to arrive.

[0303] 6.5.2 Memory Constraint Subband Data Structure

[0304] This data-structure is required to efficiently store subbandcoefficients, in memory, during the rendering algorithm. This isrequired since the coefficients are represented in long integer(lossless coding mode) or floating-point (lossy coding mode) precisionwhich typically require more memory than pixel representation (1 byte).In lossy mode, the coefficients at the client side 110 are representedusing floating-point representation, even if they were computed at theserver side 120 using an integer implementation. This will minimizeround-off errors.

[0305] At the beginning of the rendering algorithm, coefficient andpixel memory strips are initialized. dyadicWidth (ROI) may be denoted asthe width of the projection of the ROI onto the resolutiondyadicResolution (ROI). For each component and resolution 1<j≦dyadic Resolution (ROI), four subband strips are allocated for the four types ofsubband coefficients: hl,lh,hh and HalBit. The coefficient strips areallocated with dimensions$\left\lbrack {{2^{j - {{dyadicResolution}{({ROI})}} - 1} \times {{dyadicWidth}({ROI})}},{{\frac{3}{2} \times {tileLength}} + \frac{maxFilterSize}{2}}} \right\rbrack$

[0306] For each component and resolution 1≦j<dyadicResolution a pixelstrip is allocated with dimensions$\left\lbrack {{2^{j - {{dyadicResolution}{({ROI})}}} \times {{dyadicWidth}({ROI})}},{{tileLength} + \frac{maxFilterSize}{2}}} \right\rbrack$

[0307] Beginning with the lowest resolution 1, the algorithm proceedswith a recursive multiresolution march from the top of the ROI to bottom(y direction. Referring to FIGS. 10 and 11, in step 1101, themultiresolution strips are filled with sub-tiles of coefficients 1050decoded from the database or read from the memory cache. From thecoefficients we obtain multiresolution pixels 1051 using an inversesubband transform step 1102 (shown in further detail in FIG. 10). Eachtime a tile of pixels at resolutions j<dyadicResolution (ROI) isreconstructed, it is written into the pixel strip at the resolution j.Each time a tile of pixels at the highest resolution dyadicResolution(ROI) is reconstructed, it is fed into the inverse color transform andresizing steps 1103,1104.

[0308] 6.5.3 Step 1101: Decoding and Memory Caching

[0309] The subband coefficients data structure described previously insection 6.5.2 is filled on a tile basis. Each such subband tile isobtained by decoding the corresponding data blocks stored in thedatabase or reading from the memory cache. The memory cache is used tostore coefficients in a simple encoded format. The motivation is this:the decoding algorithm described previously in section 5.2 iscomputationally intensive and thus should be avoided whenever possible.To this end the rendering module uses a memory cache 111 where subbandcoefficient are stored in very simple encoded format which decodes veryfast. For each required subband tile, the following extraction procedureis performed, described in FIG. 25, beginning at step 2501. In step2502, if no data blocks are available in the database for the subbandtile, its coefficients are set to zero (step 2503). In step 2504, if thetile's memory cache storage is updated, namely it stores coefficients inthe same precision as in the database, then the coefficients can beefficiently read from there (step 2505). In step 2506, the lastpossibility is that the database holds the subband tile in higherprecision. Then, the tile is decoded down to the lowest available bitplane using the algorithm previously described in section 5.2 and thecached representation is replaced with a representation of this higherprecision information.

[0310] 6.5.4 Step 1102: Inverse Lossless Wavelet Transform

[0311] This is an inverse step to step 603 performed in the server (see7.1.5). Following FIG. 21 we see that four “extended” subbandcoefficient sub-tiles of length tileLength/2+maxFilterSize at theresolution i are read from the coefficient strips data structure andtransformed to a tile of pixels at the next higher resolution usinglosslessWaveletdTransformType(i). If i+1<dyadicResolution (ROI), thetile of pixels obtained by this step is inserted into the pixel memorystrip at the resolution i+1. If i+1=dyadicResolution (ROI) the tile ofpixels is processed by the next step of color transform.

[0312] Remark: Recall from 5.1.1 that the “half bits” are initialized aszeros, therefore the inverse step is well defined even if their “real”value is not available in the client yet.

[0313] 6.5.5 Step 1103: Inverse Color Transform

[0314] This is an inverse step to step 603 performed at the server 120.It is performed only for tiles of pixels at the resolutionhighestSubbandResolution. At this stage, all of the pixels of each suchtile are in the outputColorSpace and so need to be transformed into adisplaying or printing color space. For example, if the original imageat the server 120 is a color image in the color space RGB, the pixelsobtained by the previous step of inverse subband transform are in thecompression color space YUV. To convert back from YUV to RGB, we use theinverse step described in FIG. 13. If the ROI is at an exact dyadicresolution, then the tile of pixels in the rendering color space iswritten into the off-screen buffer. Else it is resized in the next step.

[0315] 6.5.6 Step 1104: Image Resize

[0316] In case the resolution of the ROI is not an exact dyadicresolution, the image obtained by the previous step must be re-sized tothis resolution. This can be accomplished using operating system imagingfunctionality. In most cases the operating system's implementation issub-sampling which produces in many cases an aliasing effect which isvisually not pleasing. To provide higher visual quality, the imagingsystem of the present invention may use the method of linearinterpolation, for example described in J. Proakis and D. Manolakis,“Digital signal processing”, Prentice Hall, 1996. The output of theinterpolation is written to the off-screen buffer. From there it isdisplayed on the screen using system device dependant methods.

[0317] 6.5.7 Step 1105: Mapping to 8-bit Screen

[0318] When luminanceMap is active mapping to 8-bit screen is performedusing the mapping function described in 6.1.1.

[0319] 7. Server Worflow

[0320] With reference to FIG. 5, the operation of the server computer120 (FIG. 1) will now be described. Initially, an uncompressed digitalimage is stored in, for example, storage 122 of the server computer 120.This uncompressed digital image may be a two-dimensional image, storedwith a selected resolution and depth. For example, in the medical field,the uncompressed digital image may be contained in a DICOM file.

[0321] Once the client computer 110 requests to view or print a certainimage, the server performs the preprocessing step 501. This step is acomputation done on data read from the original digital image. Theresults of the computation are stored in the server cache device 121.After this fast computation a “ready to serve” message is sent from theserver to the client containing basic information on the image.

[0322] In step 502, the server receives an encoded stream of requestsfor data blocks associated with a ROI that needs to be rendered at theclient. The server then decodes the request stream and extracts therequest list.

[0323] In step 503, the server reads from cache or encodes data blockassociated with low resolution portions of the ROI, using the cachedresult of the preprocessing stage 501.

[0324] If the ROI is a high-resolution portion of the image, the server,in step 504, reads from cache or performs a “local” and efficientversion of the preprocessing step 501. Specifically, a local portion ofthe uncompressed image, associated with the ROI, is read from thestorage 122, processed and encoded. In step 505, the data that wasencoded in steps 503-504 is progressively sent to the client in theorder it was requested.

[0325] 7.1 Step 501: Preprocessing

[0326] The preprocessing step is now described with respect to FIG. 6.The preprocessing algorithm's goal is to provide the fastest response tothe user's request to interact with the image. Once this fastcomputational step is performed, the server is able to provide efficient“pixel-on-demand” transmission of any client ROI requests that willfollow. In most cases the first ROI is a view of the full image at thehighest resolution that “fits” the viewing device. The preprocessingalgorithm begins with a request for an uncompressed image that has notbeen processed before or has been processed but the result of thisprevious computation has been deleted from the cache. As explained, thisunique algorithm replaces the possibly simpler procedure of encoding thefull image into some progressive format. This latter technique willprovide a much slower response to the user's initial request then thetechnique described below. At the end of the algorithm a “ready to serveROI of the image” message is sent to the client containing basicinformation on the image. While some of this information, imagedimensions, original color space, resolution etc., is available to theuser of the client computer, most of this information is “internal” andrequired by the client to formulate ROI request lists (§6.2) andprogressively render (§6.5). Next we describe in detail thepreprocessing algorithm.

[0327] 7.1.1 Preprocessing Parameters TABLE 7 Variable MeaninglosslessMode If true, preprocessing will prepare data that can be usedfor lossless transmission. subbandTransformType The framework allows theability to select a different subband transform for each resolution ofeach image. The technical term is: non-stationary transforms.numberOfResolutions The number of resolutions in the Multiresolutionstructure calculated for the image. Typically,${numberOfResolutions} = {{\log_{2}\left( \sqrt{ImageSize} \right)}.}$

jumpSize A number in the range [0, numberOfResolutions − 1]. Thepreprocessing stage computes only the top lower part of the image'smultiresolution pyramid of the size numberOfResolutions − jumpSize.tileLength Typically = 64. Tradeoff between time and coding performance.nTilesX (j) Number of subband tiles in the x direction at the resolutionj nTilesX (j) Number of subband tiles in the y direction at theresolution j inputColorSpace Color space of uncompressed original image.outputColorSpace Transmitted color space used as part of the encodingtechnique. numberOfComponents Number of components in OutputColorSpace.threshold (c, j) A matrix of values used in lossy compression. Thesubband coefficients of the component c at the resolution j withabsolute value below threshold (c, j) are considered as (visually)insignificant and set to zero.

[0328] Given an input image, the parameters described in Table 7 arecomputed or chosen. These parameters are also written into a header sentto the client 110 and are used during the progressive rendering step 405(see section 6.5, described previously). The important parameters toselect are:

[0329] 1. losslessMode: In this mode, progressive transmission of imagestakes place until lossless quality is obtained. Choosing this moderequires the preprocessing algorithm to use certain reversible wavelettransforms, and can slow down the algorithm.

[0330] 2. subbandTransformType(j): The (dynamic) selection of waveletbasis (as described, for example, in I. Daubechies, “Ten lectures onwavelets”, Siam, 1992) is crucial to the performance of the imagingsystem. The selection can be non-stationary, meaning a differenttransform for each resolution j. The selection is derived from thefollowing considerations:

[0331] a. Coding performance (in a rate/distortion sense): This isobviously required from any decent subband/wavelet transform.

[0332] b. Approximation of ideal low pass: It is favorable to choose atransform such that low resolutions of the image will be of high visualquality (some filters produce poor quality low resolutions even beforeany compression takes place).

[0333] c. Fast transform implementation: Can the associated fasttransform be implemented using lifting steps (as described, for example,by I. Daubechies and W. Sweldens, “Factoring wavelet transforms intolifting steps”, J. Fourier Anal. Appl., Vol. 4, No. 3, pp. 247-269,1998), using only integer shifts and additions, etc. Some good examplesare the Haar and CDF transforms (1,3), (2,2)*** described in I.Daubechies, “Ten lectures on wavelets”, Siam, 1992.

[0334] d. Fast low pass implementation: A very important parameter,since together with the parameter jumpSize, it determines almost all ofthe complexity of the algorithm. For example, the CDF (1,3) is in thisrespect the “optimal” transform with three vanishing moments. Since thedual scaling function is the simple B-spline of order 1, its low pass issimple averaging. Thus, the sequence of CDF transforms, using theB-spline of order 1 as the dual scaling function, but with wavelets withincreasing number of vanishing moments are in some sense optimal in thepresent system. They provide a framework for both real time response andgood coding efficiency.

[0335] e. Lossless mode: IflosslessMode is true we must choose thefilters from a subclass of reversible transforms (see, for example,“Wavelet transforms that map integers to integers”, A. Calderbank, I.Daubechies, W. Sweldens, B. L. Yeo, J. Fourier Anal. Appl., 1998).

[0336] f. Low system I/O: If the network 130 in FIG. 1 connectingbetween the Image residing on the storage 122 and the imaging server 120is slow, the bottleneck of the preprocessing stage (and the wholeimaging system for that fact) might be simply the reading of theoriginal image. In such a case a transform may be chosen with a lazysub-sampling low pass filter that corresponds to efficient selectivereading of the input image. Many interpolating subband transforms withincreasing number of vanishing moments can be selected to suit thisrequirement. However, this choice should be avoided whenever possible,since it conflicts with (a) and (b).

[0337] g. Image type: If the type of the image is known in advance, anappropriate transform can be chosen to increase coding efficiency. Forexample: Haar wavelet for graphic images, smoother wavelet for real-lifeimages, etc. In the graphic arts field, there are numerous cases ofdocuments composed of low resolution real-life images and highresolution graphic content. In such a case, a non-stationary sequence oftransforms may be chosen: Haar for the high resolutions and a smootherbasis starting at the highest resolution of a real-life image part. Incase of low system I/O (f), a non-stationary choice of interpolatingtransforms of different orders is required.

[0338] 3. jumpSize: This parameter controls the tradeoff between fastresponse to the user's initial request to interact with the image andresponse times to subsequent ROI requests. When jumpSize is large, theinitial response is faster, but each computation of a region of interestwith higher resolution than the jump might require processing of a largeportion of the original image.

[0339] 4. InputColorSpace: The input color spaces supported in losslessmode are:

[0340] a. Grayscale: For grayscale images

[0341] b. RGB

[0342] 5. outputColorSpace: The following are color spaces which performwell in coding:

[0343] a. Grayscale: For grayscale images

[0344] b. YUV: for viewing color images

[0345] Referring to Table 7 [LP], losslessMode is set to true. Threshold(c,j) is not in use, since in lossless mode, there is no thresholding.The rest of the variables have the same meaning as in the lossyalgorithm.

[0346] 7.1.2 Memory Constraint Multiresolution Scan Data Structure

[0347] Most wavelet coding algorithms have not addressed the problem ofmemory complexity. Usually the authors have assumed there is sufficientmemory such that the image can be transformed in memory from the timedomain to a wavelet frequency domain representation. It seems theupcoming JPEG2000 will address this issue, as did its predecessor JPEG.The preprocessing algorithm also requires performing subband transformson large images, although not always on the fall image, and thusrequires careful memory management. This means that the memory usage ofthe algorithm is not of the order of magnitude of the original image, asdescribed in J. M. Shapiro, “An embedded hierarchical image coder usingzero-trees of wavelet coefficients”, IEEE Trans. Sig. Proc., Vol. 41,No. 12, pp. 3445-3462, 1993.

[0348] Given an uncompressed image we allocate the following number ofmemory strips

numberOfComponents×(numberOfResolutions−jumpSize)

[0349] Of sizes

[2^(−(numberOfResolutions−j))imageWidth,3×tileLength/2+maxFilterSize]

[0350] For 1≦j≦numberOfResolutions−jumpSize−1 and

[image Width,tileLength+2×maxFilterSize]

[0351] For j=numberOfResolutions−jumpSize

[0352] That is, the memory usage is proportional to2^(−jumpSize)×imageWidth. Each such strip stores low-pass coefficientsin the color space outputColorSpace at various resolutions.

[0353] Referring to FIG. 6, during the preprocessing stage, theresolutions are scanned simultaneously from start to end in the ydirection. For each color component and resolution, the correspondingstrip stores low-pass coefficients or pixels at that resolution. Thecore of the preprocessing algorithm are steps 604-607, where tiles ofpixels of length tileLength+2×maxFilterSize are read from the memorystrips and handled one at a time. In step 604 the tile is transformedinto a tile of length tileLength containing two types of coefficientdata: subband coefficients and pixels at a lower resolution. The subbandcoefficients are processed in steps 605-606 and are stored in the cache.The lower resolution pixels are inserted in step 607 to a lowerresolution memory strip. Whenever such a new sub-tile of lowerresolution pixels is inserted into a strip, the algorithm's memorymanagement module performs the following check: If the part of thesub-tile exceeds the current virtual boundaries of the strip, thecorresponding first lines of the strip are considered unnecessaryanymore. Their memory is (virtually) re-allocated for the purpose ofstoring the sub-tile's new data.

[0354] 7.1.3 Step 601: Lossless Color-transform

[0355] This step is uses the conversion formula described in FIG. 13.This step must be performed before step 602, because the lossless colorconversion is non-linear. 7.1.4 Step 602: Lossless Wavelet Low Pass

[0356] The motivation for the low pass step is explained in § 6.1.4 inthe above-cited Ser. No. 09/386,264, which disclosure is incorporatedherein by reference. In a lossless mode there are a few emphasis thatare represented here.

[0357] In step 602, the low pass filters of the transformssubbandTransformType(j), numberOfResolutions−jumpSize<j≦numberOfResolutions, are used to obtain a low resolution strip at theresolution numberOfResolutions−jumpSize (as can be seen in FIG. 26).Typically, it is required to low pass about 2^(jumpSize) lines of theoriginal image 1010 to produce one line of the low resolution image. Thelow pass calculation is initiated by a read of tiles of pixels from thememory strips performed in step 604. Whenever there is an attempt toread missing low resolution lines, they are computed by low passing theoriginal image and inserted into the memory strip. The insertionover-writes lines are no longer required, such that the algorithm ismemory constrained. In the case where a non-linear color transform tookplace in the previous step 601, the results of that transform arelow-passed.

[0358] In the lossless mode of operation, the jumpSize parameter definesthe number of lossless wavelet low pass steps should be done. A singlelow pass step is the same for Haar and CDF (1,3) and defined by thefollowing two stages (taken from (3.20) and (3.22):${X\text{-}{{direction}:{s(n)}}} = \left\lfloor \frac{{x\left( {2n} \right)} + {x\left( {{2n} + 1} \right)}}{2} \right\rfloor$

[0359] Y-direction: s(n)=x(2n)+x(2n+1).

[0360]

[0361] Namely, in a 2D representation, the low pass step is defined by$\begin{matrix}{{l\quad {l\left( {m,n} \right)}} = {\left\lfloor \frac{{x\left( {{{2m} + 1},{{2n} + 1}} \right)} + {x\left( {{2m},{{2n} + 1}} \right)}}{2} \right\rfloor + {\left\lfloor \frac{{x\left( {{{2m} + 1},{2n}} \right)} + {x\left( {{2m},{2n}} \right)}}{2} \right\rfloor.}}} & (7.1)\end{matrix}$

[0362] ForjumpSize=1 and jumpSize=2 (other sizes practically are notneeded), the server performs these steps efficiently (almost like thelossy algorithm) by a single operation that simulates exactly jumpSizelow pass steps defined in (7.1). As noticed from (7.1), the simplicityof the formula makes filters such as Haar and CDF (1,3) “optimal” in therespect of low pass efficiency.

[0363] 7.1.5 Step 603: Forward Lossless Wavelet Transform

[0364] In Step 603 we perform one step of an efficient local losslesswavelet transform (§3), on a tile of pixels at the resolution1≦j≦numberOfResolutions−jumpSize. The type of transform is determined bythe parameter lossless WaveletTransformType(j). As described in FIG. 1,the transform is performed on an “extended” tile of pixels is of lengthtileLength+2×maxFilterSize (unless we are at the boundaries), readdirectly from a multi-resolution strip at the resolution j+1. The outputof the step is a lossless subband tile composed of wavelet coefficientsincluding Half bit coefficients and low resolution coefficients/pixels.The transform step is efficiently implemented in integers as describedin §3.

[0365] The subband transform of step 603 outputs three types of data:scaling function (low pass), wavelet (high pass) and Halfbits. Thewavelet coefficients are treated in step 604 while the scaling functioncoefficients are treated in step 605.

[0366] Remark: Tiles of pixels which are located on the boundariessometimes need to be padded by extra rows and/or columns of pixels, suchthat they will formulate a “full” tile of length tileLength.

[0367] 7.1.6 Step 604: Variable Length Encoding and Storage

[0368] In step 604, the subband coefficients that are calculated in step603 are variable length encoded and stored in the cache 121. IfmaxBitPlane (tile)=0 we do not write any data. Else we loop on thecoefficient groups {coef(2×i+x,2×j+y)}_(x,y=0,1). For each such group wefirst write the group's variable length length (i,j) usinglog₂(maxBitPlane(tile)) bits. Then for each coefficient in the group wewrite length (i,j)+1 bits representing the coefficient's value. Theleast significant bit represents the coefficient's sign: if it is 1 thenthe variable length encoded coefficient is assumed to be negative. TheHalfBit subband coefficients are written in one-bit per coefficient.

[0369] 7.1.7 Step 605: Copying Low Pass Coefficients into theMultiresolution Strip Structure

[0370] In step 503, unless losslessMode is true, the subbandcoefficients calculated in step 604 are quantized. This procedure isperformed at this time for the following reason: It is required that thecoefficients computed in the previous step will be stored in the cache121. To avoid writing huge amounts of data to the cache, somecompression is required. Thus, the quantization step serves as apreparation step for the next variable length encoding step. It isimportant to point out that the quantization step has no effect oncompression results. Namely, the quantization step is synchronized withthe encoding algorithm such that the results of the encoding algorithmof quantized and non-quantized coefficients are identical.

[0371] A tile of an image component c at the resolution j is quantizedusing the given threshold threshold (c,j): for each coefficients x, thequantized value is └x/threshold(c,j)┘. It is advantageous to choose theparameters threshold (c,j) to be dyadic such that the quantization canbe implemented using integer shifts. The quantization procedureperformed on a subband tile is as follows:

[0372] 1. Initialize maxBitPlane(tile)=0.

[0373] 2. Loop over each group of four coefficients{coef(2×i+x,2×j+y)}_(x,y=)0,1. For each such group initialize a variablelength parameter length (i,j)=0.

[0374] 3. Quantize each coefficient in the group coef (2×i+x,2×j+y)using the appropriate threshold.

[0375] 4. For each coefficient, update length (i,j) by the bit plane bof the coefficient, where the bit plane is defined by

|coef(2×i+x,2×j+y)|∈[2^(b)threshold(c,j),2^(b+1)threshold(c,j))

[0376] 5. After processing the group of four coefficients, use the finalvalue of length (i,j) to update maxBitPlane(tile) by

maxBitPlane(tile)=max (maxBitPlane(tile),length (i,j))

[0377] 6. At the end of the quantization step, store the valuemaxBitPlane(tile) in the cache 121.

[0378] Note that for subband tiles located at the boundaries we can setto zero subband coefficients that are not associated with the actualimage, but only with a padded portion. To do this we take into accountthe amount of padding and the parameter maxFilterSize. The motivationfor the “removal” of these coefficients is coding efficiency.

[0379] 7.2 Step 502: Decoding the Request Stream

[0380] This is the inverse step of section 6.3. Once the request streamarrives at the server 120, it is decoded back to a data block requestlist. Each data structure the type representing a group of requesteddata blocks is converted to the sub-list of these data blocks.

[0381] 7.3 Step 503: Encoding Low Resolution Part of ROI

[0382] Step 503 is described in FIG. 7. It is only performed wheneverthe data blocks associated with low-resolution subband tile are notavailable in the server cache 121.

[0383] Step 701 is the inverse step of step 604 described in §7.1.6. Inthe preprocessing algorithm subband tiles of lower resolution, that isresolutions lower than numberOfResolutions−jumpSize, were stored in thecache using a variable length type algorithm. For such a tile we firstneed to decode the variable length representation. The algorithm usesthe stored value maxBitPlane(tile).

[0384] 1. If maxBitPlane(tile)=0, then all the coefficients are set tozero including the HalfBit subband.

[0385] 2. If maxBitPlane(tile)=1, then all the coefficients are set tozero, and the HalfBit subband coefficient are read bit by bit fromcache.

[0386]3. Else, as performed in 2, the HalfBit subband coefficient areread bit by bit from cache, and we perform the following simple decodingalgorithm:

[0387] For each group of four coefficients{coef(2×i+x,2×j+y)}_(x,y=)0,1, we read log₂(maxBitPlane(tile)) bitsrepresenting the variable length of the group.

[0388] Assume the variable length is length (i,j). For each of the fourcoefficients we then read length (i,j)+1 bits. The least significant bitrepresents the sign. The reconstructed coefficient takes the value:$\left( {{readBits}\operatorname{>>}1} \right) \times \left\{ \begin{matrix}{- 1} & {{{{readBits}\&}1} = 1} \\1 & {{{{readBits}\&}1} = 0}\end{matrix} \right.$

[0389] In step 702 we use the encoding algorithm described in §5.1 toencode the requested data blocks associated with the extracted subbandtile.

[0390] 7.4 Step 504: Processing High Resolution Part of ROI

[0391] Step 504 is described in FIG. 8. In case we have used jumpSize>0in step 501 and the resolution of the ROI>numberOfResolutions−jumpSize,we are sometimes required to perform a local variation of thepreprocessing step described in §7.1. Whenever the server receives arequest list of data blocks we check the following. If a data block hasbeen previously computed (present in the cache 121) or is associatedwith a low resolution subband tile data block then it is either simplyread from the cache or handled in step 503. Else, the coordinates of thedata block are used to find the “minimal” portion of the ROI that needsto be processed. Then, a local version of the preprocessing algorithm isperformed for this local portion. The difference here is that step 804replaces Variable Length coding step 604 of the preprocessing algorithmby the encoding algorithm given in §5.1.

[0392] 7.5 Step 505: Progressive Transmission of ROI

[0393] In the final step, the encoded data tiles are sent from theserver 120 to the client 110, in the order they were requested. In manycases, data blocks will be empty. For example, for a region of theoriginal image with a constant pixel value, all of the correspondingdata blocks will be empty, except for the first one that will containonly one byte with the value zero representing maxBitPlane(tile)=0. Fora low activity region of the image, only the last data blocksrepresenting higher accuracy will contain any data. Therefore, to avoidthe extra side information, rectangles of empty data blocks arecollected and reported to the client 110 under the restriction that theyare reported in the order in which they were requested. For blockscontaining actual data, only the data block's size in bytes need bereported, since the client 110 already knows which data blocks to expectand in which order.

[0394] The present invention has been described in only a fewembodiments, and with respect to only a few applications (e.g.,commercial printing and medical imaging). Those of ordinary skill in theart will recognize that the teachings of the present invention may beused in a variety of other applications where images are to betransmitted over a communication media

What is claimed:
 1. A system for transmitting a digital image over acommunication network, comprising: (a) an image storage device forstoring a digital image; (b) a client computer coupled to thecommunication network, wherein the client computer generates andtransmits across the communication network coordinates defining a regionof interest within the digital image; (c) a server computer, coupled tothe communication network and the image storage device, wherein theserver computer performs the steps of: (i) pre-processing the digitalimage through a lossless wavelet transformation; (ii) receiving thecoordinates from the client computer; and (iii) progressivelytransmitting to the client computer the region of interest within thedigital image defined by the coordinates.
 2. The system of claim 1,wherein the server computer progressively transmits the region ofinterest to a selected quality threshold.
 3. The system of claim 1,wherein the server computer progressively transmits the region ofinterest to lossless quality.
 4. The system of claim 1, wherein theclient computer reverse transforms the region of interest received fromthe server computer, to form a lossless reproduction of the digitalimage.
 5. The system of claim 4, wherein client computer displays thelossless reproduction of the digital image on a web browser resident onthe client computer.
 6. The system of claim 1, wherein the servercomputer performs the pre-processing step through a lossless wavelettransformation comprising two non-identical one-dimensional transforms.7. A system for transmitting a digital image over a communicationnetwork, comprising: (a) an image storage device for storing a digitalimage; (b) a client computer coupled to the communication network,wherein the client computer generates and transmits across thecommunication network coordinates defining a region of interest withinthe digital image; (c) a server computer, coupled to the communicationnetwork and the image storage device, wherein the server computerperforms the steps of: (i) pre-processing the digital image through alossless wavelet transformation; (ii) generating wavelet coefficientscorresponding to the digital image; (iii) generating half-bit flags,each half-bit flag corresponding to a wavelet coefficient; (iv)receiving the coordinates from the client computer; (v) progressivelytransmitting to the client computer the region of interest within thedigital image defined by the coordinates, based on the waveletcoefficients and the corresponding half-bit flags.
 8. The system ofclaim 7, wherein the server computer progressively transmits the regionof interest to a selected quality threshold.
 9. The system of claim 7,wherein the server computer progressively transmits the region ofinterest to lossless quality.
 10. The system of claim 7, wherein theclient computer reverse transforms the region of interest received fromthe server computer, to form a lossless reproduction of the digitalimage.
 11. The system of claim 10, wherein client computer displays thelossless reproduction of the digital image on a web browser resident onthe client computer.
 12. The system of claim 7, wherein the servercomputer performs the pre-processing step through a lossless wavelettransformation comprising two non-identical one-dimensional transforms.13. A method for transmitting a digital image from a server computer toa client computer, wherein the client computer generates coordinatesdefining a region of interest within the digital image, and transmitsthe coordinates to the server computer, the method comprising the stepsof: (a) storing a digital image within an image storage device; (b)pre-processing the digital image through a lossless wavelettransformation; (c) receiving the coordinates at the server computerfrom the client computer; and (d) progressively transmitting from theserver computer to the client computer the region of interest within thedigital image defined by the coordinates.
 14. The method of claim 13,wherein the server computer progressively transmits the region ofinterest to a selected quality threshold.
 15. The method of claim 13,wherein the server computer progressively transmits the region ofinterest to lossless quality.
 16. The method of claim 13, furthercomprising the step of reverse transforming at the client computer theregion of interest received from the server computer, to form a losslessreproduction of the digital image.
 17. The method of claim 16, furthercomprising the step of displaying at the client computer the losslessreproduction of the digital image on a web browser resident on theclient computer.
 18. The method of claim 13, wherein the pre-processingstep is performed through a lossless wavelet transformation comprisingtwo non-identical one-dimensional transforms.
 19. A method fortransmitting a digital image from a server computer to a clientcomputer, wherein the client computer generates coordinates defining aregion of interest within the digital image, and transmits thecoordinates to the server computer, the method comprising the steps of:(a) storing a digital image within an image storage device; (b)pre-processing the digital image through a lossless wavelettransformation; (c) receiving the coordinates at the server computerfrom the client computer; and (d) generating wavelet coefficientscorresponding to the digital image; (e) generating half-bit flags, eachhalf-bit flag corresponding to a wavelet coefficient; (f) receiving thecoordinates at the server computer from the client computer; and (g)progressively transmitting from the server computer to the clientcomputer the region of interest within the digital image defined by thecoordinates, based on the wavelet coefficients and the correspondinghalf-bit flags.
 20. The method of claim 19, wherein the server computerprogressively transmits the region of interest to a selected qualitythreshold.
 21. The method of claim 19, wherein the server computerprogressively transmits the region of interest to lossless quality. 22.The method of claim 19, further comprising the step of reversetransforming at the client computer the region of interest received fromthe server computer, to form a lossless reproduction of the digitalimage.
 23. The method of claim 22, further comprising the step ofdisplaying at the client computer the lossless reproduction of thedigital image on a web browser resident on the client computer.
 24. Themethod of claim 19, wherein the pre-processing step is performed througha lossless wavelet transformation comprising two non-identicalone-dimensional transforms.
 25. A server computer for transmitting adigital image to a client computer, wherein the client computergenerates and transmits to the server computer coordinates defining aregion of interest within the digital image, the server computercomprising: (a) an image storage device for storing the digital image;(b) a processor for performing the steps of: (i) pre-processing thedigital image through a lossless wavelet transformation; (ii) receivingthe coordinates from the client computer; and (iii) progressivelytransmitting to the client computer the region of interest within thedigital image defined by the coordinates.
 26. The system of claim 25,wherein the server computer progressively transmits the region ofinterest to a selected quality threshold.
 27. The system of claim 25,wherein the server computer progressively transmits the region ofinterest to lossless quality.
 28. The system of claim 25, wherein theclient computer reverse transforms the region of interest received fromthe server computer, to form a lossless reproduction of the digitalimage.
 29. The system of claim 28, wherein client computer displays thelossless reproduction of the digital image on a web browser resident onthe client computer.
 30. The system of claim 25, wherein the servercomputer performs the pre-processing step through a lossless wavelettransformation comprising two non-identical one-dimensional transforms.31. A server computer for transmitting a digital image to a clientcomputer, wherein the client computer generates and transmits to theserver computer coordinates defining a region of interest within thedigital image, the server computer comprising: (a) an image storagedevice for storing the digital image; (b) a processor for performing thesteps of: (i) pre-processing the digital image through a losslesswavelet transformation; (ii) generating wavelet coefficientscorresponding to the digital image; (iii) generating half-bit flags,each half-bit flag corresponding to a wavelet coefficient; (iv)receiving the coordinates from the client computer; (v) progressivelytransmitting to the client computer the region of interest within thedigital image defined by the coordinates, based on the waveletcoefficients and the corresponding half-bit flags.
 32. The system ofclaim 31, wherein the server computer progressively transmits the regionof interest to a selected quality threshold.
 33. The system of claim 31,wherein the server computer progressively transmits the region ofinterest to lossless quality.
 34. The system of claim 31, wherein theclient computer reverse transforms the region of interest received fromthe server computer, to form a lossless reproduction of the digitalimage.
 35. The system of claim 34, wherein client computer displays thelossless reproduction of the digital image on a web browser resident onthe client computer.
 36. The system of claim 31, wherein the servercomputer performs the pre-processing step through a lossless wavelettransformation comprising two non-identical one-dimensional transforms.